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Abstract 

The  reduction  of  an  aircraft’s  radar  cross  section  can  increase  its  survivability  in 
hostile  airspace  by  making  it  more  difficult  to  locate  and  track  by  enemy  radar.  Replacing 
articulated  flight  control  surfaces  with  adaptive  controls  will  reduce  surface  discontinuities, 
and  enhance  low  observability.  Actuation  of  the  aerodynamic  surfaces  is  achieved  by  an 
electric  field  applied  to  PZT  actuators  embedded  in  the  top  and  bottom  skins,  creating 
differential  strain  and  shear  in  the  host  substrate.  This  creates  torsion  about  the  elastic 
axis,  and  a  change  in  the  wing  lift  coefficient.  The  torsion  of  the  designed  baseline  UAV’s 
wing  torquebox  was  modeled  in  the  presence  of  a  full  complement  of  air-loads  by  extending 
the  Bredt-Batho  theorem.  This  was  accomplished  through  modifying  Libove’s  method, 
using  a  thin-walled,  linearly  elastic,  fully  anisotropic,  trapezoid  cross-section  beam.  The 
linear  tip  twist  angles  due  to  a  uniform  cross-sectional  moment  were  verified  using  the 
isotropic  Bredt-Batho  theorem,  and  published  anisotropic  results  by  applying  isotropic, 
then  anisotropic  laminate  elastic  properties.  The  isotropic  solutions  were  within  3.1%; 
the  anisotropic  results  were  within  6.9-10.9%  of  the  published  angles.  The  PZT  actuation 
of  the  host  structure  was  achieved  by  substituting  the  PZT-composite  laminate  elastic 
properties  into  the  derived  solution  and  inducing  strain  and  shear  of  the  PZT  lamina  by 
applying  an  electric  field,  without  the  presence  of  external  forces  or  moments.  Using  two 
different  PZT  laminae,  the  angular  twist  as  a  function  of  the  host  lamina  orientation  angle 
and  applied  voltage  was  recorded.  The  amount  of  twist  ranged  between  0.03-0.39  degrees, 
and  0.12-1.04  degrees  for  the  AFC  and  G-1195  PZT  laminae  respectively. 
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MODELING  PIEZOCERAMIC  TWIST 

ACTUATION 

IN  SINGLE-CELL  ANISOTROPIC 
TORQUE  BOX 

OF  LOW-OBSERVABLE  UAV  WING 


I.  Introduction 

Mathematics  up  to  the  present  day  have  been  quite  useless  to  us  in  regard  to  flying. 


From  the  14th  Annual  Report  of  the  Aeronautical  Society  of  Great  Britain,  1879 


1.1  Background 

1.1.1  Aerial  Reconnaissance.  Strategic  and  tactical  reconnaissance  have 
been  the  sources  of  intelligence  information  for  military  commanders  since  the  first  hu¬ 
man  armed  conflict.  One  form  of  today’s  military  intelligence  gathering  is  overflight  by 
reconnaissance  vehicles,  such  as  a  satellite  or  aircraft.  These  sources  provide  detailed 
and  accurate  information  that  is  mostly  independent  of  local  weather  (Synthetic  Aperture 
Radar,  or  infrared  photography);  however,  they  cannot  always  meet  real-time  demand. 
One  might  need  to  wait  hours  before  the  satellite  is  in  position  or  until  the  long-range  air¬ 
craft  from  a  far  away  base  is  launched.  Today’s  battle  commanders  and  combat  controllers 
demand  up  to  the  minute  information  on  enemy  positions  as  well  as  real-time  battle  dam¬ 
age  assessment  in  almost  all  weather  conditions.  The  tactical  unmanned  aerial  vehicles 
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(UAVs)  can  provide  those  pieces  of  information.  They  are  deployed  close  to  the  conflict 
areas,  can  reach  their  targets  quickly,  loiter  above  the  target  area  for  extended  periods, 
and  they  can  also  be  quickly  turned  around  for  new  missions. 

The  UAVs  for  strategic  reconnaissance  in  our  current  inventory  are  divided  into 
three  categories.  The  conventional,  medium  altitude  endurance  UAV,  Predator  (Tier  II) 
is  designed  to  provide  24-hour,  near  continuous,  on-station  surveillance  with  a  500  nm 
operational  radius  using  simultaneous  carriage  of  electro-optical  (EO),  infrared  (IR)  and 
synthetic  aperture  radar  (SAR)  sensors,  at  an  altitude  of  25,000  feet.  Demonstrated  system 
capability  can  provide  20  hours  total  flight  time  at  13,000  ft  [7]. 


Figure  1.1  MAE  Predator  UAV 


The  conventional,  high-altitude  endurance  (HAE)  UAV,  Global  Hawk  (Tier  11+)  is 
designed  to  provide  24  hour,  on-station  surveillance  with  a  3000  nm  radius  using  EO,  IR 
and  SAR  sensors  at  an  altitude  of  over  50,000  ft.  The  higher  altitude  and  longer  operational 
radius  allows  greater  survivability  and  operational  flexibility  [7]. 
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Figure  1.2  HAE  Global  Hawk  on  the  dry  lakebed  of  Edwards  AFB,  CA. 

The  low-observable  (LO),  high-altitude  endurance  UAV,  Darkstar  (Tier  III)  was 
designed  for  low  observability  and  optimized  for  moderate  endurance,  high-threat  recon¬ 
naissance  missions  in  which  coverage  is  more  important  than  range  or  endurance.  It  is 
equipped  with  either  EO  or  SAR  sensors  and  flies  at  altitudes  greater  than  45,000  ft,  with 
on-station  flight  time  of  eight  hours  and  a  500  nm  operational  radius  [7].  This  program, 
however,  was  cancelled  in  1998. 

There  are  two  tactical  reconnaissance  UAVs  in  the  current  Army  and  Navy  inventory. 
The  Hunter  is  a  short-range  UAV  operated  by  the  Army,  and  is  designed  for  a  maximum 
altitude  of  15,000  ft  with  a  maximum  range  of  144  nm.  The  on-station  endurance  is  11.6 
hours.  The  Pioneer  is  a  short-range,  ship-launched  tactical  UAV  operated  by  the  Navy.  It 
flies  at  a  maximum  altitude  of  15,000  ft,  to  a  maximum  combat  radius  of  100  nm  and  has 
an  on-station  endurance  of  five  hours  [7]. 
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Figure  1.3  LO  DarkStar  UAV  in  flight 


Only  Global  Hawk  uses  jet  propulsion  (Darkstar  was  cancelled  in  1998),  and  none 
employs  low-observable  (stealth)  technology.  This  factor  was  a  major  contributor  to  los¬ 
ing  four  Predators  (three  lost  to  surface-to-air  missiles  /SAMs/),  six  Hunters  (four  lost  to 
SAMs),  and  four  Pioneers  (three  lost  to  SAMs)  during  Operation  Allied  Force  [7].  While 
fortunately  all  of  these  losses  were  calculated  in  dollars  and  not  in  human  lives,  the  vulner¬ 
ability  of  low-flying,  non-st ealthy  UAVs  remains  a  concern.  This  is  why  the  Air  Combat 
Command  (ACC),  the  Air  Force  Research  Laboratory  (AFRL)  and  the  Defense  Advanced 
Research  Projects  Agency  (DARPA)  are  currently  conducting  a  joint  unmanned  combat 
aerial  vehicle  (UCAV)  Advanced  Technology  Demonstration  (ATD)  program.  In  Phase  II 
of  the  program  (Engineering  and  Manufacturing  Development,  or  EMD  of  the  acquisition 
process)  Boeing  will  develop  two,  low-observable  technology  (stealth)  8,000  lb,  tailless,  34 
ft  wingspan,  jet-powered  UCAVs  [7]. 


1.1.2  Issues  in  Low- Observability.  We  have  seen  that  most  of  the 
tactical  UAVs  in  the  Kosovo  and  Bosnia  operations  were  lost  due  to  enemy  fire.  We  also 
know  that  an  aircraft’s  survivability  in  a  hostile  environment  can  be  greatly  enhanced 
by  employing  stealth  technology,  which  reduces  the  vehicle’s  radar  cross  section  (RCS) 
compared  to  that  of  a  similar  category  and  size  of  aircraft.  There  are  a  number  of  ways 
to  reduce  an  aircraft’s  RCS,  among  them  are  the  use  of  radar  wave  absorbing  materials 
and  paints,  radar  wave  directing  surfaces,  enclosed  ordinance,  stores  or  equipment,  and 
hidden  surface  discontinuities  (such  as  weapons,  engines,  antennas).  In  order  to  reduce 
the  RCS  of  the  aircraft’s  wing,  the  structural  discontinuities  caused  by  the  aerodynamic 
control  surfaces  (ailerons,  flaps,  slats  or  spoilers)  need  to  be  reduced  or  eliminated.  The 
flight  control  of  the  aircraft  then  will  have  to  be  accomplished  by  adaptive  control  using 
adaptive  aerodynamic  control  surfaces.  These  surfaces  (wings,  horizontal  and  vertical 
tail)  will  have  to  be  actuated  so  that  they  can  generate  the  necessary  increase  in  the  lift 
coefficient  required  for  a  particular  type  of  maneuver. 

1.1.3  Adaptive  Controls.  The  idea  of  using  aeroelastic  control  by  deform¬ 
ing  lifting  surfaces  was  first  implemented  by  the  Wright  brothers.  They  achieved  control 
of  the  Wright  flyer  by  warping  the  end  of  the  lifting  surfaces  by  the  means  of  cables  and 
pulleys.  This  control  scheme,  however,  proved  awkward  due  to  the  complexity  of  its  de¬ 
sign,  and  was  ineffective  at  higher  dynamic  pressures  due  to  the  relatively  primitive  nature 
of  the  aircraft  itself.  Today,  with  more  sophisticated  tools  and  materials  available,  it  is 
possible  to  realize  the  benefits  of  adaptive  control,  particularly  through  the  use  of  strain 
actuation  [16].  An  adaptive  airfoil  equipped  with  strain  actuators  can  then  be  used  for 
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control  by  inducing  twist,  or  camber  in  the  wing,  rather  than  using  articulated  control 
surfaces  [16].  The  induced  twist  and  camber  can  then  be  regulated  to  produce  the  desired 
aerodynamic  forces  and  moments. 

1.2  Induced  Strain  Actuation 

Induced  strain  actuation  is  the  process  by  which  actuation  strain  in  the  elements  of  a 
structure  induces  deformation  of  the  overall  structure.  One  of  the  most  commonly  used  in¬ 
duced  strain  actuation  methods  is  piezoelectricity,  in  which  an  applied  electric  field  creates 
strain  in  the  piezoelectric  material.  Shape-memory  metal  alloys  (SMA),  torque-tubes,  mag- 
netostrictives,  electrostrictives,  and  piezoelectrics,  in  particular  Lead-Zirconium-Titanate 
(PZT)  piezoceramics  have  been  used  as  actuators.  It  has  many  advantages  over  other  types 
of  actuation  because  the  actuators  are  easily  integrated  within  the  load  bearing  structures 
by  either  bonding  it  onto  the  surface  or  embedding  it  in  the  structural  element.  [10]. 

A  number  of  studies  have  been  published  on  the  analysis  of  strain  actuation  using 
PZTs  as  actuators.  The  majority  of  the  research  involved  the  deflection  control  of  rectan¬ 
gular  plates  and  beams  using  surface  bonded  piezoceramic  elements.  Crawley  and  Lazarus 
(1989)  have  developed  analysis  techniques  for  strain-actuated,  plate-like  adaptive  struc¬ 
tures,  showing  that  induced  strain  actuation  is  an  effective  means  for  controlling  those 
structures.  The  strain  actuation  was  achieved  by  using  composite  materials,  piezoceram¬ 
ics,  and  shape  memory  alloys.  Their  wind  tunnel  experiments  demonstrated  that  sufficient 
static  aeroelastic  control  can  be  developed  using  the  adaptive  structures  [15].  Batra  (1995) 
illustrated  the  use  of  PZTs  as  sensors  and  actuators  to  control  the  deflection  of  the  cen- 
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troid  of  a  rectangular  plate  subjected  to  a  uniformly  distributed  load.  To  analyze  the 
problem  they  developed  a  finite  element  model  employing  four-noded  Lagrangian  elements 
[4].  Crawley  and  Anderson  (1990)  modeled  induced  strain  actuation  of  beam-like  compo¬ 
nents  of  intelligent  structures.  They  constructed  a  Bernoulli-Euler  beam  model  of  surface 
bonded,  and  embedded  actuators,  that  included  the  extension  and  bending  of  the  actuator. 
In  order  to  resolve  the  uncertainty  whether  the  surface  bonded  actuator  they  developed 
agrees  with  the  Bernoulli-Euler  beam  model,  they  constructed  a  finite  element  model  that 
included  the  extension,  bending  and  shear  in  the  actuator  structure  [8]. 

Because  closed  form  solutions  cannot  be  found  for  the  majority  of  induced  strain 
actuation  problems,  most  of  the  work  done  with  actual  aerodynamic  control  surfaces  dealt 
with  the  problem  experimentally  or  used  approximate  solution  methods  [10].  Forster  and 
Yang  (1998)  investigated  the  effects  of  using  piezoceramic  actuation  in  the  flutter  control 
of  wing  boxes.  Assuming  that  the  piezoceramic  actuation  can  achieve  the  desired  change 
in  the  natural  frequency  of  the  wing  box,  they  demonstrated  that  significant  savings  in 
wing  structural  weight  can  be  achieved  [12],  Kudva,  et.  al.  (1996)  investigated  the  static 
adaptive  control  of  aerodynamic  surfaces  with  embedded  shape  memory  alloys  (SMA) 
providing  antagonistic  actuation.  Their  tests  at  NASA  Langley  in  1996  demonstrated 
an  8%  increase  in  lift  due  to  a  wing  twist  of  1.25  deg  [14].  Romeo,  et.  al.  investigated 
the  linear  and  nonlinear  angle  of  twist  of  rectangular  composite  laminate  torqueboxes, 
and  contrasted  the  predicted  values  of  twist  with  experimental  results.  They  found  good 
correlation  between  the  theoretical  analysis  and  the  experimental  results  by  considering 
the  effects  of  the  non-linear  shear  modulus  caused  by  incomplete  diagonal  shear  stress  [24]. 
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Relatively  few  papers  approached  the  problem  using  direct  analytical  methods  with¬ 
out  resorting  to  finite  element  methods.  The  wing  structure  —  composed  of  spars  and  skin 
—  is  generally  modeled  as  a  closed-,  single-,  or  multiple-cell,  thin- walled  beam  section. 
The  torsion  of  the  wing  can  be  analyzed  as  torsion  of  the  closed  section  thin- walled  beam. 
Badir  (1995)  developed  a  variationally  and  asymptotically  consistent  theory  in  order  to 
predict  the  static  response  of  anisotropic,  thin-walled,  two-cell  beams  subjected  to  exten- 
sional  loads,  torsion  and  bending.  One  of  the  major  advantages  of  his  study  was  that  the 
displacement  field  was  not  assumed  a  priori,  and  it  emerges  as  the  result  of  the  asymptotic 
analysis  of  the  shell  energy  [2].  Libove  (1988)  generalized  the  Bredt-Batho  formula  for 
doubly-symmetric,  anisotropic  shells.  In  order  to  reduce  the  complexity  of  the  solution, 
Libove  simplified  the  problem  to  circular,  and  doubly  symmetric,  high  slenderness  ratio, 
bi-convex  tubes  of  unit  thickness. 

1.3  Scope  and  Approach 

This  thesis  investigates  the  possibility  of  eliminating  the  articulated  ailerons  by  equip¬ 
ping  the  aircraft  with  an  adaptive  wing  that  uses  piezoceramic  torsion  actuation  for  roll 
control.  For  this  effort,  a  specific  UAV  design  using  conceptual  aircraft  design  meth¬ 
ods  is  developed  in  order  to  obtain  the  wing  structural  parameters,  such  as  spar  heights, 
torque-box  areas  and  skin  thickness.  The  overall  design  is  used  to  build  a  scale  model  for 
experimental  scattering  and  RCS  testing  using  the  AFRL/SNA  RCS  evaluation  facility. 
See  reference  [25]  for  details  on  the  theoretical  RCS  computer  simulation  model  and  its 
comparisons  with  the  experimental  results.  The  selected  wing  design  was  a  thin-walled, 
closed,  three-section,  box-beam  construction.  The  shear  flow  solution  for  the  torsion  of 
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a  single-cell,  isotropic  thin-walled  box  beam  is  going  to  be  developed  using  the  Bredt- 
Batho  method.  The  technique  is  going  to  be  expanded  for  the  torsion  solution  of  a  closed, 
three-cell,  isotropic  box-beam.  The  classical  theory  will  then  be  generalized  to  a  single-cell 
isotropic  and  anisotropic  model  in  order  to  account  for  elastic  couplings  present  in  modern 
laminated  composite  wing  designs.  To  validate  the  accuracy  of  the  generalized  anisotropic 
solution,  the  generalized  model  isotropic  result  is  compared  with  the  Bredt-Batho  results. 
The  generalized  anisotropic  result  will  be  compared  to  the  results  obtained  by  Romeo. 
To  account  for  the  PZT  actuator  lamina  incorporated  in  the  structure,  a  continuous  PZT 
element  was  assumed  to  be  embedded  in  the  top  and  bottom  surfaces  of  the  single-cell 
box  beam.  The  material  elastic  properties  of  the  PZT-composite  element  are  compared  to 
those  of  the  composite  alone  construction  skin.  The  angular  tip  displacement  due  to  tor¬ 
sion  is  compared  to  that  of  the  composite  only  construction  beam.  The  structure  is  strain 
actuated  by  applying  an  electric  field  to  the  PZT  actuator  lamina,  and  the  generated  twist 
angles  will  be  compared  to  the  angles  generated  by  uniform  torsion. 

1-4  Overview 

The  current  presented  background  information  on  current  issues  in  UAV  and  low- 
observable  technology.  It  also  provided  a  review  on  some  of  the  analytic  and  experimental 
work  done  on  the  subject  of  bending  and  torsion  of  beams,  plates,  and  aerodynamic  surfaces 
using  piezoceramic  strain  actuation.  Chapter  II  will  detail  the  conceptual  design  of  the 
proposed  low-observable  UAV  in  order  to  obtain  the  wing  structural  parameters,  so  that 
the  wing  structural  geometry  is  based  on  a  specific  and  realistic  UAV  design.  It  will 
identify  the  wing  torquebox  geometric  dimensions  for  one  specific  gross  weight  UAV,  as 


1-9 


well  as  those  of  the  simplified  torquebox  in  order  to  aid  in  the  calculations  presented  in 
Chapters  III  and  IV.  This  chapter  will  also  include  several  other  UAV  designs  based  on 
maximum  gross  takeoff  weight  and  wing  aspect  ratios  in  order  to  facilitate  possible  future 
extensions  of  this  study  to  include  trade  studies  of  takeoff  weight  and  piezoelectric  strain 
actuation. 

Chapter  III  includes  the  shear  flow  solution  for  the  torsion  of  a  closed,  single-cell, 
thin-walled,  isotropic  box-beam,  and  the  extension  of  Bredt-Batho  theorem  to  a  closed- 
section,  three-cell,  isotropic  beam.  The  required  applied  moment,  the  resulting  shear  flows 
in  the  individual  structural  members,  and  the  tip  angle  of  twist  is  detailed. 

Chapter  IV  derives  the  generalization  of  the  Bredt-Batho  theorem  to  a  fully  aniso¬ 
tropic,  single-cell,  asymmetric,  thin- walled  box  beam  through  the  modification  of  Libove’s 
solution  for  the  torsional  displacement  of  closed-section,  doubly-symmetric,  single-cell,  cir¬ 
cular  tubes.  To  predict  the  accuracy  of  the  solution,  the  method  is  applied  to  the  torquebox 
geometry,  designed  in  Chapter  II,  by  considering  it  as  an  isotropic  single-cell  beam,  and  the 
solution  obtained  is  then  compared  to  that  of  the  isotropic  thin- walled  beam  of  Chapter  III. 
The  accuracy  of  the  anisotropic  torsion  model  is  compared  to  the  values  published  in  ref¬ 
erence  [24]  for  the  torsion  of  an  anisotropic,  rectangular  box  beam  subjected  to  uniform 
cross-sectional  moments. 

Chapter  V  extends  the  results  of  Chapter  IV  to  a  fully  anisotropic,  single-cell,  thin- 
walled  box  beam  with  PZT  lamina  embedded  in  the  center  of  the  top  and  the  bottom 
skins.  The  material  properties  of  the  PZT-composite  will  modify  the  constitutive  relations 
of  Chapter  IV,  and  the  resulting  tip  twist  angle  due  to  uniform  cross-sectional  moment 
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acting  on  the  single-cell  box  beam  will  be  compared  to  those  calculated  from  Chapter  IV. 


Piezoelectric  strain  is  then  applied  to  the  PZT-composite  construction,  anisotropic  torque- 
box,  and  the  resulting  tip  torsional  displacement  is  recorded  as  a  function  of  applied  voltage 
and  substrate  lamina  angle. 

Chapter  VI  will  detail  the  results  obtained  from  Chapters  III  through  V.  Chapter  VII 
will  conclude  and  summarize  the  work  included  in  this  thesis,  and  will  make  recommenda¬ 
tions  as  to  the  direction  of  further  studies  or  theses  on  this  subject. 
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II.  Conceptual  Design  of  Tactical  UAV 

I  began  to  realize  that  there  might  be  something  after  all  to  Newton’s  Laws. 

Robert  H. Goddard,  1902 


2.1  Operational  Requirements 

The  operational  requirements  for  the  proposed  low-observable  UAV  were  established 
in  order  to  provide  a  design  that  bridges  the  gap  between  the  high-endurance,  high- weight, 
long-range  Global  Hawk,  and  the  short-range,  light-weight  Predator.  To  establish  realistic 
wing  and  torque  box  dimensions,  some  operational  requirements  were  set,  and  six  baseline 
UAVs  were  designed:  three  different  takeoff  weights,  each  with  two  different  aspect  ratios. 
All  other  design  parameters  (such  as  taper  ratio,  payload  to  gross  weight  percentage,  etc.) 
were  held  constant.  The  aircraft  operational  requirements  were  set  as: 

•  Gross  Takeoff  Weight  of  5,000,  10,000,  and  15,000  lb 

•  Useful  Payload  of  500,  1,000,  and  1,500  lb 

•  Single-Engine  Jet  Propulsion 

•  Wing  Aspect  Ratio  of  10  and  12 

•  Mission  Altitude  of  25,000  ft 

•  Endurance  of  10  hours  or  more 

•  Range  of  2,000  NM,  or  more 

•  Use  of  High  Endurance  Airfoil  NLF(1)0215 
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2.2  Conceptual  Design  Philosophy 

During  conceptual  design  the  engineer  uses  approximations  to  initially  estimate  the 
aircraft  structural  parameters  and  performance.  For  the  initial  estimates,  conceptual  de¬ 
sign  equations  are  used,  which  in  turn  are  iterated  until  the  required  design  parameters 
are  satisfied.  More  iterations  can  improve  the  design,  but  they  also  considerably  increase 
the  time  and  level  of  effort  required  to  achieve  the  required  performance  and  structural 
parameters.  The  purpose  of  the  designs  for  this  study  was  to  obtain  structural  as  well 
as  performance  parameters  that  are  realistic  to  the  given  applications  at  hand,  and  not 
necessarily  to  design  an  aircraft  that  is  ready  to  fly.  In  the  case  of  the  design  of  a  complete 
aircraft,  the  conceptual  design  is  followed  by  the  phases  of  preliminary,  then  detailed  de¬ 
sign.  These  phases  further  refine  the  performance  and  structural  parameters  by  considering 
the  effects  of  subsystems  (cargo,  engine  size,  landing  gear  fitting,  etc.)  on  the  overall  design 
of  the  aircraft,  and  by  considering  the  fulfillment  of  civilian  (Federal  Aviation  Regulations 
or  FARs  set  by  the  Federal  Aviation  Administration)  or  Military  Standards  (MILSTDs). 
For  more  detailed  information  about  conceptual,  preliminary  and  detailed  aircraft  design, 
the  interested  reader  is  referred  to  reference  [22]. 

2.3  Initial  Aircraft  Sizing 

The  initial  sizing  of  the  low-observable,  tailless,  single-engine,  jet  UAV  was  proposed 
in  order  to  facilitate  the  manufacturing  of  a  scaled  model  to  be  used  as  part  of  comparative 
RCS  study,  as  well  as  part  of  studying  actuated  wing  torsion.  The  initial  sizing  was  carried 
out  in  order  to  provide  approximate  aircraft  structural  data  that  is  suitable  for  the  purposes 
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Table  2.1  I 

Design  Parameters 

Range  (nm) 

2,000 

Wing  Dihedral  Angle  (deg) 

2 

Endurance  (hr) 

10 

Specific  Fuel  Consumption  (SFC) 

0.4 

Aspect  Ratio  (AR) 

12 

Payload  (Wpl)(lb) 

500 

Taper  Ratio  (A) 

0.40 

Fuel  Weight  Fraction  ( Wr /) 

0.45 

Wing  Thickness  Ratio  (t) 

0.15 

Empty  Weight  Fraction  (Wre) 

0.45 

Cruise  Speed  (kts) 

180 

Pre-Flight  Fraction  (fr{) 

0.97 

Stall  Speed  (kts) 

80 

Climb  Weight  Fraction  ( fr 2) 

0.985 

Oswald  Efficiency  (e) 

0.8 

Cruise  Weight  Fraction  ( fr 3) 

0.870 

of  these  studies,  and  it  was  not  needed  to  be  further  refined  by  the  methods  of  preliminary 
and  detailed  aircraft  design.  Therefore,  only  two  sets  of  iterations  were  performed  resulting 
in  sufficiently  accurate  numbers  for  the  purpose  of  the  study. 


2.3.1  Design  Parameters.  In  order  to  achieve  the  desired  performance, 
several  conceptual  structural  ratios  and  performance  factors  were  assumed  based  on  ac¬ 
cepted  initial  structural  sizing  ratios.  Two  aspect  ratios  (AR)  of  10  and  12  were  set,  and 
an  Oswald  efficiency  factor  (which  accounts  for  the  span  efficiency  due  to  the  non-elliptical 
lift  distribution,  as  well  as  for  the  variations  of  the  parasite  drag  with  lift)  of  0.80  was  cho¬ 
sen  for  the  low-speed,  high  L/D,  low-wing  aircraft  configuration.  The  fuel-to-gross  weight 
(FGW)  ratio  of  0.45  was  assumed  on  the  basis  of  FGW  ratio  of  0.35  of  light  business  jets 
over  12,000  lb  gross  takeoff  weight  (GTW).  The  empty-to-gross  weight  (EGW)  ratio  of  0.45 
was  assumed  on  the  basis  of  extensive  composite  structure,  and  the  EGW  of  0.55  of  all- 
metal  construction,  light  business  jets  over  12,000  lb  GTW.  The  specific  fuel  consumption 
(SFC)  of  0.4  was  assumed  on  the  basis  of  older,  less  efficient  engines  SFC  of  0.5-0. 6.  Fur¬ 
ther  refinement  and  optimization  of  these  parameters  via  preliminary  and  detailed  design 
is  suggested  as  topic  of  future  research.  The  design  parameters  are  shown  in  Table  (2.1). 
Using  the  design  parameters  and  the  aircraft  weight  fractions  from  Table  (2.1),  the  various 
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aircraft  and  fuel  weights  can  be  calculated: 


Wq  =  Gross  Weight  =  - —  - — — —  (2.1) 

9  y  1  -Wre-Wrf 

We  =  Empty  Weight  =  WreWg  (2.2) 

Wf  =  Fuel  Weight  =  WrfWg  (2.3) 

Wi  =  Preflight  Weight  =  fr\Wg  (2.4) 

W2  =  Climb  Weight  =  fr2Wx  (2.5) 

Wz  =  Cruise  Weight  —  frzW2  (2.6) 


2.3.2  Aerodynamic  Data.  The  initial  design  calculated  the  wing  surface 
area  required  for  the  given  cruise  conditions  in  standard  atmosphere  by  accounting  for  the 
fuselage  form  (lift)  factor  F,  the  decrease  of  lift-curve  slope  due  to  compressibility  (Mach) 
effects,  finite  wing  and  wing  sweep  effects,  and  the  reduction  of  wing  area  by  the  fuselage. 
The  wing  area  for  the  given  stall  speed  was  calculated  and  was  found  to  be  greater  than 


that  required  for  the  cruise  conditions.  Because  high-lift  devices  such  as  leading  edge 


Table  2.2  Aerodynamic  Data  for  N 


Climax 

1.7 

ao  (1/rad) 

2i r 

o-stall  (deg) 

13 

ao L  (deg) 

-5.0 

aii  (deg) 

0 

Cd 

0.008 

/F-0215  Airfoil 


slats  or  slotted  flaps  are  not  considered  due  to  the  low  RCS  considerations,  the  wing  area 
necessary  for  the  slow  flight  (stall  condition)  was  selected  for  the  design.  Although  wing 
twist  (wash-out)  improves  on  the  stall  characteristics  of  an  aircraft  with  tapered  wings, 
this  design  did  no  select  twist  (as  a  function  of  span)  in  order  to  simplify  the  modeling  of 
the  RCS,  ease  of  manufacturing,  and  simplified  torsion.  The  aerodynamic  data  obtained 
from  published  airfoil  data  for  the  high-altitude,  long  endurance  NLF(1)0215  airfoil  are 
listed  in  Table  (2.2),  where  Cimax  is  the  maximum  profile  (2-D)  lift  coefficient,  astaii  is  the 
profile  (2-D)  stall  angle  of  attack  (AOA),  a*  is  the  wing  incidence  angle  (arbitrarily  set), 
ao  is  the  profile  (2-D)  lift-curve  slope,  aoz,  is  the  profile  zero-lift  AOA,  and  Cd  is  the  2-D 
profile  drag  coefficient  at  ao l-  Using  these  parameters,  the  approximate  finite  wing  (3-D) 
lift-curve  slope  aw ,  the  Mach  number  correction  factor  /3,  and  the  lift  curve  slope  corrected 
for  compressibility  effects  77  can  be  calculated  via  Eqns.  (2.7)  through  (2.9). 


&W 


ARao 


A  P  _l  2AR-\- 4 
+  AR+2 


(2.7) 


(2.8) 


Tj  = 


ao 


180 


O  7T 

ZP 


(2.9) 
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Table  2.3  Structural  Reference  Values 


Reference  Wing  Span,  bref  (ft) 

42 

Reference  Wing  Sweep,  Q,^ref)  (rad) 

4.08  •  10"3 

Reference  Wing  Area,  5re/  (ft2) 

151 

Exposed  Wing  Area  Ratio,  Ratio 

0.825 

Fuselage  Width,  Wfuse  (ft) 

7  bref 

Maximum  Wing  Lift  Coefficient,  C^max 

1.52 

2.3.3  Reference  Parameters.  In  order  to  be  able  to  continue  with  some 
of  the  performance  parameter  calculations,  it  was  necessary  to  set  some  of  the  structural 
and  performance  parameters  as  reference  parameters.  These  values  were  then  updated 
through  repeated  iterations  until  the  desired  convergence  was  obtained.  The  reference 
values  necessary  for  the  further  steps  in  the  design  are  listed  in  Table  (2.3).  These  reference 
values  were  used  in  calculating  the  fuselage  form  (lift)  factor  F,  the  lift  curve  slope  for  the 
entire  aircraft  aw,  the  aircraft  cruise  lift  coefficient  Cicrmse,  the  cruise  wing  loading  Lw, 
and  the  total  wing  surface  required  for  cruise,  Scruise.  These  were  calculated  using 


F  =  1.07 


(2.10) 


2nAR 


a  in  — 


4  + 


(  AB?P 

VI3- 


)( 


1  + 


(tan(Qr 

W 


+  2 


&exp  p  ft 

Sref  180 


(2.11) 


Chemise  —  Q"w  {^cruise  -  Qol) 


(2.12) 


_  2  r1 

2  '  P  ’  vcr  *  ^ Lcruise 


(2.13) 


2-6 


Scruise  =  (2-14) 

■L'W 

where  p  is  the  air  density  at  mission  altitudes,  and  Sexp  is  the  exposed  wing  surface  area 
not  covered  by  the  fuselage,  given  by 

$ * (>xp  —  i^idtzo  *  Sj 'sf  (2.15) 

The  wing  area  for  cruise  ( Scruise )  obtained  from  Eqn.  (2.14)  will  be  compared  to  the  wing 
area  calculated  from  stall  conditions. 

2.3.4  Wing  Area.  Now  that  some  of  the  structural  as  well  as  some  of 
the  aerodynamic  performance  parameters  were  determined  for  cruise  conditions  in  Sec¬ 
tion  2.3.3,  the  same  were  calculated  for  stall  speed  at  standard  sea  level  conditions.  The 
values  obtained  in  this  present  section  were  then  substituted  into  the  reference  values  of 
Section  2.3.3,  and  were  iterated.  For  the  stall  conditions  at  sea  level,  the  wing  loading, 
wing  surface,  wing  span,  wing  chord  lengths  as  well  as  thickness  at  the  root  and  the  tip, 
wing  sweep  angle,  and  the  maximum  lift  coefficient  due  to  the  wing  sweep  were  calculated: 

Lyj  ~  g  ’  PO  ^ cr  '  ^ Lmax  (2.16) 


b  —  yj AR  ■  S stall 


(2.18) 


Croot  — 


25 

fo(l-t-A) 


Ctip  —  A  '  Croot 


(2.19) 


t root  —  ^  ’  Croot  % 


tip 


i  *  Cfip 


(2.20) 


ft  =  atari 


2  •  Croot  \ 

b  ) 


(2.21) 


CL  max  —  0-9  Climax  *  COs(A) 


(2.22) 


where  t  is  the  wing  profile  thickness  from  Table  (2.1),  the  factor  0.9  in  Eqn.  (2.22)  is  the 
generally  accepted  scaling  factor,  and  A  is  the  quarter-chord  sweep  angle  given  by 


A  =  atari 


3 


(1  A)  croot 
2b 


(2.23) 


Using  the  values  obtained  from  Eqns.  (2.16)  through  (2.23),  we  can  define  the  reference 
values  for  the  exposed  wing  area  Sexp  and  the  parameter  Ratio ,  discussed  in  Section  2.3.3. 
These  were  calculated  by 


x  — 


— tan(fl) 
14 


(2.24) 
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cfuse 


—  Croot  2# 


(2.25) 


Srov  —  2 


2(cfn 


+  Cr0ot ) 


wfv 


(2.26) 


Sexp  =  s-  Scov  (2.27) 

Ratio  =  (2.28) 

O 

where  x  is  the  decrease  in  the  chord  length  due  to  the  effects  of  the  wing  taper,  Cfuse  is 
the  wing  chord  length  at  the  fuselage-wing  connection,  SCOv  is  the  wing  area  covered  by 
the  fuselage,  Sexp  is  the  exposed  wing  area  not  covered  by  the  fuselage,  and  finally  the 
parameter  called  Ratio  is  the  value  used  in  the  initial  calculations  in  Section  2.3.3. 

2 . 3. 5  Drag  and  Power  Required .  In  order  to  calculate  the  endurance  and 
range  of  the  aircraft,  the  drag  and  power  required  at  various  flight  phases  and  conditions 
must  be  determined.  The  results  can  also  be  used  during  the  preliminary,  then  detailed 
design  to  narrow  the  selection  of  suitable  engines,  and  later  to  select  the  appropriate  engine 
for  the  aircraft  as  well  as  for  the  mission.  To  find  the  drag  estimate  during  cruise  and  stall, 
we  calculate  the  lift  coefficients  for  both  conditions.  The  drag  coefficient  estimate  at  cruise 
and  at  stall  will  be  the  sum  of  the  parasitic  and  induced  drag  coefficients.  The  drag 
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coefficient  at  cruise  is  given  by 


C Demise 


c2 

r 1  i  ^ Lcruise 

Cd  +  ~^AR 


(2.29) 


where  e  is  the  Oswald  efficiency  factor  (see  Section  2.3.1)  and  AR  is  the  wing  aspect  ratio. 
The  cruise  lift  coefficient  in  Eqn.  (2.29)  is  given  by 


r  _  2  W2 

^ Lcruise  —  o  n 

Pvcr  S 

Then,  the  drag  estimate  at  cruise  is  obtained  from 


(2.30) 


D 


cruise 


cruise 


s 


(2.31) 


In  unaccelerated,  straight  and  level  flight,  the  drag  is  equal  to  the  thrust  required.  The 
power  required  is  calculated  by  multiplying  the  drag  (thrust  required)  by  the  cruise  velocity. 


Pr  —  D cruise^  cr  —  PR^cr 


(2.32) 


The  drag,  thrust  and  power  required  for  stall  conditions  can  be  readily  calculated  using 
Eqns.  (2.29)  through  (2.32).  For  stall  conditions,  the  values  of  Ci^max ,  and  Cpcruise  are 
substituted,  wherever  appropriate. 

2.3.6  Range  and  Endurance.  For  a  jet  powered  aircraft,  the  specific  fuel 
consumption  (SFC)  is  defined  as  the  weight  of  fuel  consumed  per  unit  thrust  per  unit  time 
[13],  Because  the  fuel  consumption  of  a  jet  airplane  depends  on  the  thrust  produced  by 


2-10 


the  engine,  the  SFC  for  a  jet  airplane  is  defined  as  the  thrust  specific  fuel  consumption 
(TSFC),  denoted  by  ct.  With  the  assumption  that  ct  is  constant  in  time,  for  the  endurance 
of  the  jet-powered  aircraft  we  get  : 


^(g)U) 


(2.33) 


Naturally,  the  lift  and  drag  coefficients  are  for  cruise  conditions,  Wg,  and  W3  are  as  defined 
in  Section  2.3.1.  The  division  by  3600  was  used  to  obtain  the  endurance  in  the  customary 
units  of  hours,  rather  than  seconds.  For  the  background  and  the  derivation  of  Eqn.  (2.33), 
the  interested  reader  is  referred  to  reference  [13]. 

The  range  of  the  jet  powered  aircraft  is  calculated  from  the  modified  Breguet  range 
formula  (see  reference  [13]).  Assuming  constant  ct,  Cl,  Cd,  and  p  at  cruise  speed  and 
conditions  we  have 


(2.34) 


Naturally,  the  lift  and  drag  coefficients  are  for  cruise  conditions,  Wg ,  and  W3  are  as  defined 
in  Section  2.3.1.  The  division  by  6076  was  used  to  obtain  the  range  in  the  customary  units 
of  nautical  miles  rather  than  units  of  feet. 


2.3.7  Initial  Sizing  Summary.  The  initial  sizing  did  not  account  for  the 
loss  of  lift  and  increased  drag  due  to  the  downforce  generated  by  the  longitudinal  control 
surfaces  (horizontal  tail,  and  trailing  edge  elevator);  the  additional  drag  generated  by  the 
fuselage,  or  the  lift  lost  due  to  positive  dihedral  and  wing  in-flight  bending;  the  change 
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Table  2.4  Initial  Sizing  Results 


Aircraft  Parameter 

Value 

Aircraft  Parameter 

Value 

Wg  (lb) 

5000 

6(ft) 

42.66 

We  (lb) 

2250 

croot  (ft) 

5.08 

Wf  (lb) 

2250 

ctip  (ft) 

2.03 

Wi  (lb) 

4850 

n  (deg) 

13.39 

W2  (lb) 

4780 

CLmaxsweep 

1.52 

wing  (1/deg) 

0.09796 

A(deg) 

6.12 

P 

0.95 

x  (ft) 

0.73 

V 

0.99 

cfuse 

3.63 

F 

1.3976 

Scov  (ft2) 

26.53 

aw  a/c  (1/deg) 

0.1108 

Sexp  (ft2) 

125.16 

CLaruise 

0.77 

Ratio 

0.82509 

Lwcruise  (Jb/ ft  ) 

38.19 

Cj)cruise 

0.0216 

S cruise  ( ft  ) 

130.93 

Demise  (lb) 

161.09 

Sexp  (ft2) 

125.16 

Pr 

48900 

Lwstall  (lb/ ft2) 

32.96 

E  (hr) 

13.73 

Sstall  (lb/ ft2) 

151.69 

R  (nm) 

2410 

in  the  maximum  lift  coefficient  due  to  the  increase  in  the  Reynolds  number  (Re).  The 
initial  sizing  neither  calculated,  nor  attempted  to  achieve  any  stability  derivatives.  As  it 
was  pointed  out  earlier,  preliminary  and  detailed  design  must  take  these  variables  into 
consideration.  Not  accounting  for  these  factors  in  the  conceptual  design  phase  does  not 
reduce  the  validity  of  the  model  intended  to  be  used  in  the  proposed  studies. 


The  results  of  the  initial  sizing  after  two  sets  of  iterations  are  now  summarized  in 
Table  (2.4).  These  values  are  presented  for  the  5,000  lb  gross  takeoff  weight  air  vehicle. 


2.4  Wing  Sizing 

2.4-1  Air  Loads  Estimates.  The  modeling  of  the  system  of  aerodynamic 
loads  acting  on  the  aircraft  at  any  given  flight  phase  and  atmospheric  condition  becomes 
increasingly  complicated  as  our  demand  for  accuracy  grows.  To  develop  a  true  represen- 
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tation  of  the  acting  forces  and  moments,  extensive  finite  element  packages  (NASTRAN, 
ASTROS,  PATRAN)  are  commercially  available.  However,  for  conceptual  design  purposes, 
calculating  the  forces  and  moments  even  to  sliderule  accuracy  was  more  than  adequate. 
As  it  was  stressed  before,  this  would  not  be  satisfactory  if  prototyping  is  imminent.  Pre¬ 
liminary  and  detailed  design  must  address  the  loads  more  accurately.  For  initial  structural 
sizing,  however,  the  aerodynamic  loads  (forces  and  moments)  acting  on  the  wing  can  be 
closely  approximated  by  the  sum  of  distributed  forces  acting  as  point  loads  and  concen¬ 
trated  moments.  Therefore,  the  maximum  lift  force  (vertical  shear)  is  closely  approximated 
by  the  expression 


L  =  (FS)(LLF)Wg  (2.35) 

where  the  variable  FS  is  the  Factor  of  Safety,  and  LLF  is  the  Limit  Load  Factor.  The  FS  is 
commonly  set  to  FS= 1.5  for  manned,  and  FS=  1.2  for  unmanned  aerial  vehicles.  The  LLF 
(also  abbreviated  by  n)  is  commonly  set  to  n= 3.8  for  normal  category,  n=4.4  for  utility, 
n— 6.0  for  aerobatic,  and  n= 2.5-3  for  transport  category  aircraft.  Fighter  aircrafts  draw 
their  own  mission  specialized  LLF  requirements.  For  definitions  on  aircraft  categories, 
please  refer  to  reference  [22].  The  maximum  bending  moment  acting  on  the  wing  will  be 
conceptually  modeled  as 


(2.36) 


where  L  is  obtained  from  Eqn.  (2.35),  and  b  is  the  wingtip  to  wingtip  span.  It  is  important 
to  note  while  sizing  the  spars  and  the  skin,  that  one  wing  carries  only  one  half  of  L 
calculated  from  Eqn.  (2.35). 

2.4.2  Spar  Web  Sizing.  In  order  to  properly  size  the  main  and  the  rear  spar 
of  the  wing,  a  location  of  25%,  and  70%  were  chosen  respectively.  The  selection  of  the  25% 
chord  location  of  the  main  spar  followed  the  generally  accepted  practice,  since  the  pitching 
moment  about  the  quarter  chord  of  the  airfoil  is  independent  of  the  angle  of  attack.  The 
70%  chord  location  of  the  rear  spar  was  chosen  as  to  lessen  the  section’s  torsional  stiffness 
by  locating  the  two  spars  as  close  as  practically  acceptable.  It  was  initially  assumed  that 
the  main  spar  and  the  rear  spar  will  carry  2/3  and  1/3  of  the  shear  load  respectively.  A 
factor  of  safety  (FS)  of  1.2  was  designated  following  the  general  practice  for  unmanned,  low 
maneuverability  vehicles.  A  limit  load  factor  (LLF)  of  2.5  was  set,  in  the  effort  of  designing 
the  structure  for  a  maximum  of  2.5  g  load  factor.  Finally,  it  was  assumed  that  the  main 
spar  and  the  rear  spar  should  carry  2/3  and  1/3  of  the  bending  moment,  respectively. 

As  a  general  procedure,  the  structural  components  were  sized  and  then  evaluated  for 
over-,  or  underdesign  by  evaluating  the  margin  of  safety  (MS).  The  MS  is  given  by 


MS  = 


CriticalLoad 

ActualLoad 


(2.37) 


where  the  critical  load  is  the  appropriate  maximum  shear,  compressive  or  tensile  force, 
or  bending  moment  allowed  to  be  acting  on  the  member  (driven  by  material  properties). 
The  actual  load  is  the  corresponding  actual  maximum  load  designed  to  be  acting  on  the 
structure. 
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In  aircraft  wing  spar  analysis  it  is  commonly  assumed  that  the  spar  caps  absorb  all  of 
the  bending  stress,  and  the  web  (that  is  extended  to  the  full  depth  of  the  spar)  carries  all 
of  the  shear  [22].  The  actual  load  situation  is  (as  always)  somewhere  in  between;  however, 
for  conceptual  design  purposes,  these  assumptions  are  more  than  sufficiently  accurate.  It  is 
also  assumed  that  the  shear  is  constant  within  the  web,  and  therefore,  the  maximum  shear 
stress  equals  the  average  shear  stress.  The  thickness  of  the  main  spar  web  was  determined 
using  reference  [21].  Using  the  maximum  design  load  L\ ,  we  have 


_  3  Lx 

vw&bli  a 


2  T s' 


hi 


(2.38) 


where  L\  is  the  maximum  load  on  the  main  spar,  Tsmax  is  the  maximum  shear  stress 
allowed  (material  properties),  and  h\  is  the  height  of  the  main  spar  web  (given  by  profile 
parameters).  The  spar  web  will  fail  in  buckling  long  before  the  actual  material  maximum 
shear  stress  is  reached.  In  order  to  account  for  the  critical  shear  buckling  stress  Fscr ,  the 
shear  web  buckling  parameter  K,  and  web  aspect  ratio  rweb  (a/b  in  reference)  was  selected 
from  reference  [22].  The  critical  shear  buckling  stress  is  given  by 

Fscr  =  KE  (2-39) 

where  E  is  the  material’s  Young’s  modulus.  The  actual  shear  stress  in  the  spar  web  is 
calculated  by 


Fsact  — 


Li 


h\ tmeb  i 


(2.40) 
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Table  2.5  Spar  Web  Sizing 


Parameter 

Main  Spar 

Rear  Spar 

Web  Height  (in) 

8.4334 

5.3035 

Web  Thickness  (in) 

0.072 

0.049 

Web  Aspect  Ratio,  a/b 

1.5 

1.5 

Shear  Web  Buckling  Parameter,  K 

11 

11 

Area  Moment  of  Inertia  (in4) 

Critical  Buckling  Stress  (psi) 

8338 

9765 

Actual  Stress  (psi) 

8234 

9620 

Margin  of  Safety 

0.013 

0.015 

Then,  as  discussed  above,  the  MS  of  the  web  is  given  by  substituting  Eqns.  (2.39)  and 
(2.40)  into  Eqn.  (2.37),  that  is 


MSwebl  =  -  1  (2.41) 

The  same  procedure  is  repeated  for  the  rear  spar  with  the  appropriate  values  of  the  maxi¬ 
mum  design  load  and  spar  height.  Initially,  the  web  thickness  underestimated  the  required 
thickness.  After  evaluating  the  critical  buckling  load,  a  few  iterations  are  necessary  in 
order  to  obtain  a  mandatory  positive  margin  of  safety.  In  designing  weight  critical  struc¬ 
tures  (such  as  aircraft),  a  1-5%  MS  subject  to  the  actual  application  at  hand  is  generally 
accepted.  For  the  main  spar  and  rear  spar,  the  MS  was  calculated  to  be  1.3%,  and  1.5% 
respectively.  The  calculated  values  from  Eqns.  (2.38)  through  (2.41),  including  the  values 
for  the  rear  spar  as  well,  are  shown  in  Table  (2.5). 

2.4 '3  Spar  Cap  Sizing .  The  spar  caps  are  usually  equal,  or  unequal  length 
angle  sections  bolted  (or  riveted)  onto  the  ends  of  the  spar  webs.  Their  purpose  is  to  carry 
the  bending  loads  (as  was  discussed  in  Section  2.4.2)  and  to  provide  attachment  surface 
for  securing  the  wing  skin.  The  procedure  for  selecting  the  appropriate  angle  section  was 
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to  iteratively  evaluate  the  necessary  area  moment  of  the  spar  (web  and  caps),  and  check 
whether  it  satisfies  the  required  area  moment  of  inertia  calculated  from  the  engineer’s 
theory  of  bending  (see  reference  [21]).  The  required  moment  of  inertia  for  bending  is 

Ire,  =  —V  (2.42) 

&max 

where  M  is  the  maximum  bending  moment  acting  on  the  spar,  omax  is  the  maximum 
material  stress  allowed,  and  y  is  the  maximum  distance  from  the  web  centerline  to  the  top 
of  the  web.  The  combined  area  moment  of  inertia  for  the  spar  is  given  by 


Lspar 


—  leaps  +  I web 


(2.43) 


where  the  area  moment  of  inertia  for  the  caps  (all  four)  is  calculated  using  the  parallel  axis 
theorem 


leaps  —  ^[1  angle  4”  A angle{y  Hangle )  ]  (2.44) 

The  angle  properties  must  be  iterated  (different  angle  sections  selected)  until  the  inequality 

Ireq  ^  1-spar  (2.45) 


is  satisfied. 
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Table  2.6  Spar  Cap  Sizing 


Parameter 

Main  Spar 

Rear  Spar 

Unequal  Length  Angle  Section 
Angle  Section  Area  (in2) 
Wangle  (in^) 

Uangle  (^n) 

Angle  thickness  (in) 

Radius  of  Gyration  (in) 

^caps  (in^) 

3-spar  (in  ) 

3-req  (in  ) 

3x2.5x3/16 

0.996 

0.907 

0.888 

3/16 

0.954 

47.77 

51.37 

44.97 

3x2.5x3/16 

0.996 

0.907 

0.888 

3/16 

0.954 

16.02 

16.63 

14.14 

The  same  procedure  was  repeated  for  both  spars,  and  the  angle  sections  were  iterated 
until  the  inequality  of  Eqn.  (2.45)  was  satisfied.  The  final  selection  of  caps  for  both  spars 
are  shown  in  Table  (2.6). 


2.4*4  Skin  Sizing .  The  theory  of  shear  flows  will  be  reviewed  in  Chapter  III, 
Section  3.1.1.  In  the  following  discussion,  it  will  be  assumed  the  reader  is  familiar  with  the 
theory  of  shearflows,  as  well  as  the  associated  stresses  and  stress  resultants.  Also,  reference 
[20]  is  a  good  resource  for  the  following  discussion. 

To  approximate  the  geometry  of  the  torquebox,  and  to  calculate  the  moment  of 
inertia  about  the  horizontal  axis,  the  average  skin  length  lave ,  and  average  spar  height  have 
is  calculated. 


lave  =  r1  (2-46) 

*22 

have  =  {r  (2'4?) 

h2 
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where  I21  and  I2 2  are  the  top  and  bottom  skin  lengths  of  the  torquebox  respectively,  hi 
and  /12  are  the  main  and  rear  spar  heights,  respectively.  By  summing  the  shearflow  con¬ 
tributions  all  around  the  perimeter  we  get 


L  hi 

Qij  —  Q(i-  1J-1)  +  2 j^cap  2 


(2.48) 


where  qij  is  the  shearflow  in  the  different  sections  of  the  box  beam,  and  L  is  the  moment 
arm  of  the  acting  vertical  shear  (assumed  acting  at  the  geometric  center  of  box  the  beam). 
The  shear  flow  qo  in  the  spar  web  can  also  be  easily  and  very  accurately  approximated  by 

_  T  +  2Aoqi 
^  2Aq  +  2Aq 

where  T  is  the  acting  moment  due  to  the  vertical  shear,  and  Aq  is  one  half  the  enclosed 
area  of  the  box  beam.  Because  uniform  torsion  assumes  the  cross  section  is  sufficiently 
braced  against  cross  sectional  deformation  and  warping,  the  rib  spacing  and  skin  panel 
aspect  ratio  (a/b)  must  be  accounted  for.  Using  the  guidelines  from  reference  [22],  the 
shear  buckling  parameter  if,  and  skin  panel  aspect  ratio  a/b  is  selected  to  be 


if  =  11  ;  a/b  =1.6 


(2.50) 


The  critical  shear  stress  for  buckling  of  the  flat  skin  is  again  given  by  Eqn.  (2.39),  and  is 
supplemented  by  the  critical  shear  stress  for  buckling  of  the  curved  skin  Fcra  given  by 


771  _  TP  I  TS  T?  ^skin 

r crit  —  “ scr  "r  -ft  l& 

r 


(2.51) 
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Table  2.7  Skin  Thickness  Sizing 


Parameter 

Value 

Average  Skin  Length,  lave  (in) 

28.536 

Average  Spar  Height,  have  (in) 

6.868 

Half  Area  of  Torque  Box,  Aq  (in2) 

98 

Area  Moment  of  Inertia,  I  (in4) 

98.852 

Approximate  Shearflow,  q\  (lb /in) 

271 

Radius  of  Skin  Curvature,  r  (in) 

30 

Flat  Plate  Critical  Shear  Stress,  Fscr  ( psi ) 

4118 

Curved  Plate  Critical  Shear  Stress,  Fcrit  (psi) 

5782 

Actual  Skin  Shear  Stress,  fs  (psi) 

5687 

Margin  of  Safety,  M S 

0.017 

where  K\  =  0.1  is  given  as  empirical  factor,  and  r  is  the  radius  of  curvature  for  the  skin 
panel.  The  skin  shear  is  then  calculated  by 


fs  = 


<?i 


tskin 


(2.52) 


where  q\  is  the  shearflow  calculated  from  Eqn.  (2.48),  and  tskin  is  the  initial  guess  for  the 
skin  thickness.  The  Eqns.  (2.51)  through  (2.52)  must  be  iterated  until 


Fcrit  ^  fs 


(2.53) 


The  measure  of  merit  is  again  going  to  be  the  margin  of  safety,  calculated  by 


MSskin 

Js 


(2.54) 


The  results  of  the  skin  thickness  calculations  and  iterations  for  the  5000  lb  gross  takeoff 
weight  UAV  are  presented  in  Table  (2.7). 
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2.5  Wing  Torque  Box  Identified 

2.5.1  Three- Cell,  Thin- Walled  Torque  Box .  In  order  to  facilitate 

the  shearflow  calculations  of  Chapter  III,  Section  3.1  and  Section  3.2,  all  of  the  geometric 
properties  of  the  wing  box  had  to  be  identified.  In  order  to  do  this  conveniently,  a  Matlab 
subroutine  called  ’AREA’  was  written  to  calculate  the  areas  of  the  first,  second,  and  third 
closed  sections  of  the  three-cell  beam.  The  subroutine  is  conveniently  called  by  a  driver 
program,  to  determine  the  total  skin  lengths  of  the  nose  section,  the  upper  and  lower 
center  section,  as  well  as  the  upper  and  lower  skin  lengths  of  the  trailing  edge  section.  The 
corresponding  enclosed  areas,  A\,  A<i,  and  A3  are  also  calculated.  It  is  also  necessary  to 
determine  the  actual  heights  hi,  and  h%  of  the  main  and  rear  spar  respectively  at  their 
corresponding  locations,  previously  determined  in  Section  2.4.2.  It  is  important  to  realize 
that  the  height  of  the  main  or  rear  spar  is  not  going  to  be  the  actual  dimensional  distance 
from  the  nose  multiplied  by  the  profile  percent  thickness.  That  would  lead  to  erroneous 
results.  To  avoid  this,  the  subroutine  INTERP  was  written,  which  is  used  by  the  subroutine 
AREA  in  a  two-step  iterative  linear  interpolation  scheme  in  order  to  accurately  locate  the 
positions  and  calculate  the  heights  of  the  main  and  rear  spars. 

The  driver  program  can  also  call  a  subroutine  called  PLOTPROFILE  that  will  ac¬ 
tually  draw  a  picture  representation  of  the  NLF(1)-0215F  high-altitude,  high-endurance 
airfoil.  The  airfoil  data  is  stored  in  the  subroutine  NLF0215.  The  subroutines  AREA,  IN¬ 
TERP,  PLOTPROFILE,  and  NLF0215  are  included  in  Appendix  A.  The  non-dimensional 
airfoil  shape  is  shown  in  Figure  (2.1).  A  dimensional  plot  of  the  NLF(1)-0215F  profile  with 
the  wing  root  as  well  as  at  the  wing  tip  is  included  in  Appendix  A. 
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Figure  2.1  The  NLF(1)-0215F  High-Endurance  Airfoil  Profile 

2.5.2  Simplified ,  Single- Cell,  Thin- Walled  Torque  Box.  Develop¬ 
ing  the  three-cell  closed  section  torque  box  in  Section  2.5.1  will  aid  in  the  calculations  of 
shear  flows,  torques,  and  twist  angles  in  Chapter  III.  The  development  of  a  more  general 
solution  in  Chapter  IV  would  not  have  been  possible  without  a  simplified  geometric  rep¬ 
resentation  of  the  center,  closed-section,  single-cell  torque  box  (with  area  A2).  This  is  due 
to  the  fact  that  the  mathematical  representation  of  the  curvatures  of  the  upper  and  lower 
surfaces  of  A2  are  not  practical  to  obtain  for  the  purposes  of  the  calculations  that  follow. 
Because  the  length  of  the  torque  box  ( croot  =  cr)  is  much  greater  than  its  maximum  thick¬ 
ness  ( tr  =  0.15cr),  and  the  radius  of  curvature  of  both  the  bottom  and  top  surfaces  (Z21, 
and  I22)  of  the  single-cell  center  section  are  much  greater  than  1  (r  >>  1),  it  is  reasonable 
to  approximate  these  surfaces  with  straight  lines. 

Though  the  lengths  Z21  and  Z22  are  not  equal  to  each  other  (as  it  was  determined 
from  the  subroutine  AREA),  their  difference  Z21  —  I22  is  much  less  than  the  length  of  the 
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torque  box;  therefore,  it  will  be  estimated  that 


d=^cr(c2-c  i)  (2.55) 

where  cr  is  the  chord  length  of  the  wing  root,  c\  is  the  percent  chord  location  of  the  main 
spar  (c\  =  0.25cr),  C2  is  the  percent  chord  location  of  the  rear  spar  (c2  =  0.75cr),  and 
d  is  the  calculated  half-distance  between  the  main  and  the  rear  spars.  Because  h\  ^  /12, 
the  path-length  correction  factor  T  was  found  that  will  relate  the  actual  lengths  of  the 
top  and  bottom  surfaces  to  the  horizontal  half-distance  (d)  between  the  spars.  Using  the 
Pythagorean  theorem  we  get 

T=Ji  +  {hir)2  (2'56) 

With  all  this  information  already  at  hand,  the  simplified  geometry  of  the  single-cell,  closed 
section  center  torque  box  of  the  5000  lb  gross  takeoff  weight,  AR=12  UAV  is  now  iden¬ 
tified  in  Table  (2.8).  Figure  (2.2)  shows  the  simplified  single-cell  thin-walled  torque-box 
dimensions.  The  equations  for  all  four  surfaces  (two  skins  and  two  spars)  in  terms  of  the 
path-length  coordinate  s  (that  is  using  x(s),  and  y(s))  will  be  identified  in  Chapter  IV, 
Section  4.1.1. 


2.6  Design  Summary  (Aircraft  baseline) 

After  completing  the  design  for  the  5000  lb,  AR= 12  configuration  aircraft,  the  gross 
weight  Wg  and  the  aspect  ratio  AR  of  the  aircraft  was  changed,  and  the  process  of  con- 
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Figure  2.2  Single-cell  Torquebox  Geometry 


Ta: 


ble  2.8  Simplified  Torque  Box  Dimensions 


Parameter 

Value 

Main  Spar  Height,  hi  (in) 

8.4334 

Rear  Spar  Height,  ^2)  (in) 

5.3035 

Main  Spar  Thickness,  t\  (in) 

0.072 

Rear  Spar  Thickness,  t<i  (in) 

0.048 

Top  Surface  Length  Z21  (in) 

27.7732 

Bottom  Surface  Length  I22  (in) 

29.2992 

Skin  Thickness,  ts  (in) 

0.048 

Spars  Half-Distance,  d  (in) 

13.7160 

Path  Length  Correction  Factor,  T 

1.00162 

ceptual  design  was  repeated  for  the  remaining  aircrafts  using  identical  design  parameters 
listed  in  Table  (2.1).  Six  overall  designs  were  accomplished  in  order  to  facilitate  future 
weight  versus  actuation  requirements  trade  studies.  The  baseline  aircraft  parameters  cal¬ 
culated  via  the  iterative  process  detailed  in  Section  2.3  are  listed  in  Table  (2.9).  The  wing 
structural  parameters  calculated  via  the  iterative  process  detailed  in  Section  2.4  are  listed 
in  Table  (2.10).  Five  scale  models  of  the  finished  design  were  manufactured  by  AFIT 
for  the  purpose  of  monostatic  RCS  testing  in  the  AFRL/SNA  facility.  The  models  were 
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Table  2.9 


Aircraft  Baseline  Designs 


Gross  Weight  (lb) 

5,000 

10,000 

15,000 

5,000 

10,000 

15,000 

Aspect  Ratio 

12 

12 

12 

10 

10 

10 

Empty  Weight  (lb) 

2250 

4500 

6750 

2250 

4500 

6750 

L/D 

29.66 

29.27 

28.86 

26.36 

25.92 

25.47 

Endurance  (hr) 

13.7 

13.5 

13.3 

12.18 

11.96 

11.77 

Range  (nm) 

2410 

2643 

2801 

2143 

2341 

2473 

Wingspan  (ft) 

42.66 

53.04 

59.12 

39.0 

48.48 

54.03 

Wing  Area  {ft2) 

151.69 

234.46 

291.24 

151.07 

235.05 

291.97 

Root  Chord  (ft) 

5.08 

6.31 

7.04 

5.57 

6.93 

7.72 

Tip  Chord  (ft) 

2.03 

2.53 

2.82 

2.23 

2.77 

3.09 

Root  Thickness  (ft) 

0.76 

ESI 

1.06 

0.84 

1.04 

1.16 

Tip  Thickness  (ft) 

0.30 

EH 

0.42 

0.33 

0.42 

0.46 

Table  2.10 


Wing  Structure  Baseline  Designs 


Gross  Weight  (lb) 

5,000 

10,000 

15,000 

5,000 

10,000 

15,000 

Aspect  Ratio 

12 

12 

12 

10 

10 

10 

ll(in) 

31.5955 

39.2456 

43.7859 

34.6431 

43.1018 

48.0153 

121  (in) 

27.7732 

34.4978 

38.4889 

30.4521 

37.8875 

42.2065 

122  (in) 

29.2992 

36.3933 

40.6036 

32.1253 

39.9692 

44.5256 

131  (in) 

18.6520 

23.1681 

25.8484 

20.4511 

25.4445 

28.3451 

132  (in) 

17.4096 

31.6249 

24.1267 

19.0888 

23.7497 

26.4571 

hi  (in) 

8.4334 

10.4754 

11.6873 

9.2469 

11.5047 

12.8161 

h2  (in) 

5.8151 

7.2349 

8.0597 

A\  (in2) 

91.9370 

141.8475 

176.5664 

110.5282 

171.0918 

212.3232 

A2  (in2) 

229.1889 

353.6102 

440.1609 

275.5348 

426.5133 

529.2986 

Az  (in2) 

38.5871 

59.5351 

74.1071 

46.3900 

71.8093 

89.1146 

tskin  (*n) 

0.058 

0.064 

0.047 

twebl  (i*1) 

0.098 

0.117 

0.075 

1 

tweb2  (i*1) 

mm 

0.079 

0.051 
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manufactured  with  and  without  articulated  ailerons  to  measure  and  compare  the  effects  of 
wing  discontinuities  on  RCS  scattering.  The  fuselage,  chime,  and  wing-body  interface  were 
designed  by  AFRL/SNA  using  low-observable  design  code.  The  results  of  the  test  will  be 
compared  to  the  scattering  code  developed  and  published  in  reference  [25].  A  picture  of 
one  finished  model  is  presented  in  Figure  (2.3). 


Figure  2.3  Scale  Model  of  Prototype  UAV  for  RCS  Testing 
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Ill .  Torsion  of  Box  Beams 


There  will  be  no  more  accidents! 


Gen.  Westover,  Assistant  Chief  Air  Corps,  ordering  to  stop  airmail  plane  crashes,  1934 


3.1  Single- Cell  Beam  Torsion 

3.1.1  Shear  Flows  in  Thin  Webs.  The  shear  flow  q  in  a  thin- walled 
structural  element  is  defined  as  the  product  of  the  shear  stress  os  at  the  thin  wall  centerline 
and  the  thickness  t  of  the  element: 


q  =  crst  (3.1) 

The  unit  of  shear  flow  is  force  per  unit  length.  The  value  of  q  in  a  closed,  single-cell  section 
subjected  to  torsion  alone  is  constant  along  the  section,  regardless  of  the  thickness  t  [18]. 
It  is  often  necessary  to  obtain  the  resultant  force  on  a  curved  web  in  which  the  shear  flow 
q  is  constant  along  the  length  of  the  web.  The  differential  element  ds  shown  Figure  (3.1) 
has  horizontal  and  vertical  components  dz  and  dy.  The  force  on  this  element  is  qds .  The 
horizontal  force  Fz  is 


Fz  =  [  qdz  =  qz  (3.2) 

Jo 
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Figure  3.1  Differential  Element  (Open  Section) 

where  z  is  the  horizontal  distance  between  the  ends  of  the  web.  The  vertical  force  Fy 
acting  on  the  element  is 


Fy  = 


(3.3) 


where  y  is  the  vertical  distance  between  the  ends  of  the  web.  While  Eqns.  (3.1)  and  (3.2) 
are  independent  of  the  shape  of  the  web,  the  total  torsional  moment  of  the  resultant  force 
depends  on  the  shape  of  the  web.  The  induced  moment  by  the  shear  flow  along  the  web  is 


M  = 


d^4  =  2  Aq 


(3.4) 


where  A  is  the  total  area  enclosed  by  the  web. 
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For  closed,  single-section  beams,  the  constant  shear  flow  q  around  the  circumference 
of  the  web  has  no  horizontal  or  vertical  resultant  force,  since  by  Eqns.  (3.1)  and  (3.2)  the 
horizontal  and  vertical  distances  between  the  end  points  of  the  web  are  zero.  The  resultant 
of  the  shear  flow  is  a  torque  (moment)  equal  to  the  applied  external  moment  M,  taken 
about  an  arbitrary  axis  perpendicular  to  the  cross  section.  Thus, 


M  =  J2  2(A4)g  =  2 Aq 


(3.5) 


where  A  represents  the  total  area  enclosed  by  the  web. 

The  angle  of  twist  per  unit  length  of  the  closed,  single-cell  web  is  given  by  the  Bredt- 
Batho  formula  for  torsional  displacement: 


1  -\  qAs 

2AG0  2s 


(3.6) 


where  Go  is  the  arbitrarily  selected  reference  shear  modulus  of  one  of  the  structural  mate¬ 
rials  used,  and 


t*  = 


G(s) 

Go 


(3.7) 


The  quantity  t*  is  the  modulus  weighted  thickness,  and  G(s)  is  the  shear  modulus  of  the 
web  section.  In  case  of  multi-material  construction,  G  is  expressed  as  G(s)  because  the 
value  of  G  becomes  a  function  of  the  circumferential  coordinate  s.  The  practical  application 
of  the  formula  implies  that  there  is  no  warping  constraint  on  the  cross  section,  and  it  is 


3-3 


adequately  stiffened  internally  against  distortion;  that  is,  the  shape  of  the  cross  section  is 
preserved,  neglecting  any  Poisson’s  ratio  effects  [19]. 


3.1.2  Torsion  of  the  Single- Cell  Section.  Revisiting  the  geometry  of 
the  torque-box  developed  in  Chapter  II,  we  can  apply  the  Bredt-Batho  formula  directly 
by  using  Eqn.  (3.6).  Because  the  shear  flow  q  is  constant  around  the  cross  section,  and  the 
material  properties  are  identical  in  all  sections  of  the  web,  we  can  derive  by  inspection  the 
formula  for  the  geometry  of  the  single-cell  section: 


0  = 


2AGq 


qh  ,  qh2  ,  ql  2  ,  qh  i 


,  ^ 

L  Go 


ts 


+ 


t‘2 


+ 


Gjl 

Go 


+ 


Go1 


(3.8) 


where  l\  and  h  are  the  lengths  of  the  top  and  bottom  skins  of  the  torque  box;  h\  and  h<i 
are  the  heights  of  the  main  and  rear  spars;  ts,  t\  and  t,2  are  the  thickness  of  the  skin,  main, 
and  rear  spar,  respectively;  and  Gs  and  Go  are  the  shear  moduli  of  the  skin  and  spars, 
respectively.  Using  the  known  applied  moment  M,  we  can  solve  for  q  from  Eqn.  (3.5) 


Q  = 


M 
2  A 


(3.9) 


By  factoring  out  q  and  substituting  Eqn.  (3.9)  into  Eqn.  (3.8),  we  obtain 


M  [4 Td  h2  hi 

AA2G  ts  t2  +  t\ 


(3.10) 


where  we  used  the  relation  l\  m  I2  =  2 Td,  where  T  and  d  are  as  defined  in  Chapter  II. 
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In  order  to  automate  the  calculations  of  Eqns.  (3.9)  through  (3.10),  a  Matlab  routine 
was  written,  and  is  included  in  Appendix  B.  The  routine  calculates  the  moment  required 
for,  and  the  shear  flow  generated  by  a  given  angle  of  twist.  The  code  can  be  easily  modified 
by  interested  readers  to  calculate  the  angle  of  twist  generated  by  a  given  concentrated  cross 
sectional  moment,  as  well  as  the  resulting  shear  flow. 


3.2  Multi- Cell  Beam  Torsion 

3.2.1  Torsion  of  the  Multi-Cell  Section.  The  solution  for  the  torsional 
displacement  of  the  thin- walled  beam  developed  in  Section  3.1  is  directly  applicable  to 
multi-cell  beams.  However,  sections  with  multiple  cells  are  statically  indeterminate.  In 
order  to  solve  for  the  angular  displacement,  we  must  enforce  the  condition  of  continuity  of 
rotation,  that  is,  the  angular  twist  of  all  the  cells  must  be  the  same: 

e1  =  e2  =  ...  =  en  (3.11) 


where  0,  is  given  by  the  now  familiar  formula 

ft  =  1 

1  2AiG(C  t* 


(3.12) 


The  modulus  weighted  thickness  t*  is  as  defined  by  Eqn.  (3.7);  however,  in  order  to  keep 
the  notation  as  clutter  free  as  possible,  t*  will  be  substituted  using  Eqn.  (3.7). 
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In  order  to  solve  for  the  shear  flows  qi  in  the  different  sections  of  the  multi-cell  beam, 
we  must  first  enforce  the  continuity  condition  of  Eqn.  (3.11).  Using  the  three  cells  of  the 
beam  shown  in  Figure  (3.2),  by  inspection  we  can  write 


Figure  3.2  Three-Cell  Beam  (Closed  Section) 


di  = 


l 

2-AiGo 


qih  Mgi  -  gg) 

Ga.+.  Qsl+  . 

Lg0Ti  Gq^1 


(3.13) 


e2  = 


i 


2A2Gq 


Q2hi  h2(q2  ~  qz)  92^22 


Gs  j. 

L  Go^2 


Go 


h 


M<?2  ~  9i ) 
Qsl+  , 

Go  Csl 


(3.14) 


Oz 


_1 _ 

2^43^0 


<73^31  Qz^Z2  h2(qz  -  q2) 

l&*3  ^2  J 


(3.15) 


where  the  structural  variables  are  as  defined  in  Chapter  II,  ti  are  the  skin  thickness  of 
the  corresponding  closed  section  (i  =  1,2,3),  and  tsi  is  the  main  and  rear  spar  thickness 
(i  =  1,  2)  respectively.  The  variables  <?i,  q2,  and  q 3  are  the  shear  flows  in  the  nose,  middle, 
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and  tail  sections  of  the  three-cell  beam,  respectively.  Enforcing  the  continuity  condition  of 
6\  =02,  and  solving  for  q\  in  terms  of  <72  and  <73  we  get  (after  some  heavy  algebra) 


9i 


a2 


hi  4.  hi.  -l  hi  4-  h± 
G*t 2  T  tS2  ^  G*t2  T  t,  1 


92  +  r^92 


/12  ^ 
i^3 


b  I  fci  I  -^1^1 
G*tl  'r  <sl  ‘r  A2t*2 


(3.16) 


Notice  that  by  arbitrarily  selecting  the  spar’s  shear  modulus  G  as  the  reference  shear 
modulus  Go,  we  can  use  the  simplification  Go  I  Go  =  1  and  Gs/Gq  =  G*.  By  enforcing 
the  next  continuity  condition  02  =  03,  then  by  substituting  for  q\  using  Eqn.  (3.16),  and 
solving  for  <73  in  terms  of  <72  >  we  obtain,  after  some  considerable  algebra 


93 


C1C2  -  frc,  -  (ft)'- 


os -3  +  & 


42 


ts  2  tslts2 


(3.17) 


where  the  constants  C\,  C2,  C3,  and  64  are  given  by 


hi  ^2,  ^22  h]_  ^2^2 

G*t2  ts  2  G*t2  tsl  -^3^52 


(3.18) 


/1  hi  Ax  hi 

G*t\  tsi  A-2tsl 


(3.19) 


A2  (  hi  I32  h 2  \ 
A3  \G*tz  GHZ  ts 2) 


(3.20) 


Ax  /  /21  ^2  I22  hi  \ 

A2  \G*t2  ts 2  G*t2  tsl) 


(3.21) 
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We  are  now  ready  to  solve  for  q\  in  terms  of  <72 ■  Let’s  name  the  constant  in  front  of 


q%  in  Eqn.  (3.17)  to  be  S3.  Take  Eqn.  (3.17)  and  substitute  it  into  Eqn.  (3.16)  to  get 


9i 


C4  + 


Al 

tsl 


_  hi 

ts2 


S3 


c2 


■92 


(3.22) 


Again,  let’s  name  the  constant  multiplying  q2  in  Eqn.  (3.22)  as  Si. 

In  order  to  solve  for  the  individual  shear  flows  in  the  sections,  we  use  the  extended 
version  of  Eqn.  (3.5),  that  is 


M  —  ^^(2  Aqi)  —  2A\q\  +  2A2q2  +  2A.3§3  (3.23) 

Because  we  solved  for  q\  and  qs  in  terms  of  q2 ,  we  can  substitute  Eqns.  (3.16)  and  (3.17) 
into  Eqn.  (3.23)  while  using  the  convenient  constants  Si  and  S3. 


M  =  2(AiSi  +  A2  +  A^S3)q2 

From  this  we  can  readily  obtain  the  shear  flow  q2 

M 

q2~  2(i415i+A2  +  A353) 

where  the  applied  moment  M  is  given  by 

M  =  GqJ*9  —  (GJ)eff6 


(3.24) 


(3.25) 


(3.26) 
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In  Eqn.  (3.26)  the  quantity  GqJ*  is  called  the  torsional  stiffness  of  the  section,  J*  is  the 
St.  Venant’s  constant  for  uniform  torsion  given  by 


§¥ds 


(3.27) 


Substituting  Eqn.  (3.7)  for  the  modulus  weighted  thickness  t*  we  get 


G°  §  GT7W7) ds 


(3.28) 


from  which  we  obtain  the  effective  torsional  stiffness 


GoJ*  =  (GJ)eff  = 


§  G(s)t(s)dS 


(3.29) 


In  order  to  find  the  effective  torsional  stiffness  of  the  three-section  beam,  we  first  need  to 
evaluate  the  effective  torsional  stiffness  of  each  section.  The  combined  effective  torsional 
stiffness  is  given  by  the  sum  of  the  individual  torsional  stiffnesses. 


[GJ)eff1 


Gst\  Gots 


(3.30) 


i  |  .  I22  1  fci 

Gs  £2  GotS2  Gst2  Gotsi 


(3.31) 


(GJ)efh  ~  j fe2 

Gsts  G$ts  GotS2 


(3.32) 
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Then  the  effective  torsional  stiffness  (GJ)e//  for  the  three-cell  beam  can  be  calculated  by 
using  Eqn.  (3.29) 


(GJ)eff  —  {GJ)eff1  +  ( GJ)eff2  +  {GJ)eff3  (3.33) 

We  now  have  arrived  at  the  final  solution  of  the  torsion  of  the  three-cell  beam  section. 
The  required  moment  (in  units  of  in-lb)  to  generate  a  given  angle  of  twist  6  (in  units  of 
radians),  and  the  resulting  shear  flows  (force  per  unit  length,  in  this  case  lb/ in)  can  be 
calculated  using  the  following  relations. 

M  =  ( GJ)effe  (3.34) 


M 

92  ~  2(A1S1  +  A2  +  A3S3) 


(3.35) 


qi  =  Siq2 


(3.36) 


<73  =  S3q2  (3.37) 

To  automate  the  calculations  of  Eqns.  (3.5)  through  (3.37),  a  Matlab  program  called 
SHEARFLOW3CELL  was  written  to  calculate  the  shear  flow  and  moment  required  to 
generate  a  given  tip  twist  angle  for  the  given  geometry  single-cell,  closed-section  box  beam. 
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The  code  was  extended  to  calculate  the  individual  shear  flows  in  the  three  sections  of  the 


closed  box  beam,  as  well  as  the  moment  required  to  generate  the  given  tip  twist  angle.  By 
running  the  program,  a  direct  comparison  can  be  made  between  the  effectiveness  of  the 
single-cell  and  the  three-cell  beam  in  resisting  torsional  moments.  The  Matlab  program  is 
included  in  Appendix  B. 
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IV.  Generalized  Torsion  Solution 


There  will  be  no  more  flying! 

Maj  Byron  Q.  Jones,  Eastern  Airmail  Zone  Commander,  responding  to  Westover,  1934 


4.1  Homogeneous  Isotropic  Single- Cell  Beam 

4- .1.1  Geometry  and  Loading.  While  the  shear-flow  solution  of  the  one, 
or  multi-cell  torque  box  is  useful  due  to  its  relative  simplicity,  it  only  applies  to  isotropic 
beams.  It  also  applies  to  specially  orthotropic  construction,  meaning  that  at  least  one 
of  the  axes  elastic  symmetry  is  parallel  to  the  longitudinal  axis  of  the  thin-walled  beam 
[17].  One  disadvantage  of  the  method  is  that  it  does  not  account  for  applied  loads  other 
than  the  applied  torque  acting  on  the  cross  section.  Therefore,  in  this  section,  we  extend 
the  Bredt-Batho  theorem  by  generalizing  the  method  developed  by  Libove  (see  ref  [17]) 
and  develop  an  analytical  solution  for  the  torsional  displacement  of  the  fully  isotropic, 
thin-walled,  linearly  elastic,  single-cell  box  beam  in  the  presence  of  a  full  complement  of 
air  loads.  The  theory  assumes  that  no  cross-sectional  warping  occurs  (shape  of  the  cross 
section  is  preserved),  and  that  the  longitudinal  strains  vary  linearly  over  the  cross  section 
(linearly  elastic).  In  order  to  evaluate  the  validity,  as  well  as  the  accuracy  of  the  solution, 
in  Section  6.1  we  will  check  the  isotropic  solution  against  the  shear-flow  solution  developed 
in  Chapter  III. 

Let  this  beam  segment  be  subjected  to  a  system  of  external  loads  consisting  of  forces 
and  moments  applied  at  the  end  center  of  the  cross  section  as  shown  in  Figure  (4.1).  These 
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forces  and  moment  may  be  the  result  of  applied  thrust  (longitudinal  or  vertical),  vertical 
and  horizontal  shear  due  to  lift  and  drag  respectively,  twisting  moment  due  to  external 
stores  or  strain  actuation,  as  well  as  bending  moment  due  to  lift  or  structural  weight.  In 


X 

tvy 

Figure  4.1  Forces  and  Moments  Applied  to  Closed  Section  [17] 

the  actual  cross  section  there  may  be  distributed  forces  and  moments;  however,  these  will 
be  approximated  by  concentrated  forces  and  moments  acting  on  the  section.  By  taking 
a  differential  segment  dz  of  the  cross  section,  the  applied  loads  will  be  assumed  constant 
within  dz ,  but  will  be  different  from  one  segment  to  the  next.  In  effect,  the  spanwise 
varying  loads  are  approximated  by  piecewise  uniform  forces  and  moments. 

For  simplicity,  only  a  single  cell  beam  is  considered.  The  beam  is  assumed  anisotropic, 
so  that  normal  stress  at  a  point  tends  to  generate  shear  strain  (7)  and  longitudinal  strain 
(e)  as  well.  The  cross-sectional  shear  flows  tend  to  produce  longitudinal  strain  and  shear 
strain.  This  interaction  between  normal  stress  and  shear  strains  is  called  elastic  coupling. 


4-2 


Before  we  can  start  on  the  solution,  we  have  to  establish  the  cross  sectional  geometry 
of  the  beam,  and  define  the  equations  of  the  walls  of  the  section.  Establishing  the  exact 
geometry  of  the  actual  airfoil’s  center  torque  box  is  not  an  impossible  task;  however,  it 
would  not  serve  a  very  useful  purpose  as  it  would  yield  prohibitively  lengthy  symbolic 
results  throughout  the  analysis.  Therefore,  a  generalized,  simplified  geometry  is  used  that 
approximates  the  actual  geometry  of  the  torque  box  as  a  generic  trapezoid,  defined  in 
Chapter  II. 

Let  the  x-y  coordinate  system  be  defined  with  its  origin  in  the  geometric  center  of 
the  trapezoid  cross  section.  In  addition,  let  the  s  coordinate  be  defined  as  the  clockwise 
’’path-length  coordinate”  around  the  perimeter  of  the  trapezoid,  starting  from  the  upper 
left-hand  corner  (s  —  0).  Then,  it  can  be  shown  that  for  the  four  sides  of  the  trapezoid 
torque  box  we  can  define  y(s)  and  x(s)  as 


y{s)  =  < 


h^±(s-Td)  +  \(h2  +  h1) 
2 Td  -  s  +  |/i2 
[^(3Td-s  +  h2)  +  \(h2  +  h1) 
s  —  4  Td  —  Ii2  —  \h\ 


Top  Side 
Right  Side 
Bottom  Side 
Left  Side 


(4.1) 


s 

T-d 


z(s)  =  < 


d 

3d  —  7f(s  —  /12) 


-d 


Top  Side 
Right  Side 
Bottom  Side 
Left  Side 


(4.2) 
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The  clockwise  path  coordinate  s  around  the  perimeter  is  piecewise  continuous.  The 
discontinuities  in  y(s)  at  the  respective  values  of  s  (at  the  corners  of  the  trapezoid;  see 
below)  also  serve  as  the  limits  of  integrations  for  all  the  integrals  developed  below.  For  the 
sake  of  simplicity,  these  limits  will  be  abbreviated  wherever  appropriate  by  the  following 
definitions.  Starting  from  the  upper  left-hand  corner  (5  =  0),  these  limits  are 


5i  =  2dT  52  —  2dT  +  /12 


53  =  4  dT  +  h-2  54  =  4  dT  +  h,2  +  h\ 


(4.3) 


4*1*2  Constitutive  Relations .  In  order  to  arrive  at  any  meaningful  and 
relatively  simple  solution,  the  walls  of  the  box  beam  are  assumed  to  be  thin  enough,  so 
that  they  can  be  viewed  as  membranes  in  plane  stress  (<733  =0).  In  order  to  keep  the 
solution  general,  the  method  will  be  applied  to  a  torque  box  made  of  walls  of  non-unit 
thickness.  The  solution  can  also  be  easily  modified  by  interested  readers  to  unit  thickness. 
This  can  be  achieved  by  simply  assuming  t  =  1,  thereby  eliminating  all  the  wall  thickness 
dependencies  in  the  equations.  We  can  also  assume  that  the  normal  stress  resultant  per 
unit  length  Ns  along  the  5  direction  is  negligible  [17].  This  way  the  state  of  plane  stress 
can  be  described  for  the  element  as  a  shear  flow  g,  and  a  tension  flow  N,  both  having 
units  of  force  per  unit  depth.  Figure  (4.2)  shows  the  differential  element  and  the  state  of 
stress  of  the  differential  element.  The  normal  strain  e  and  shear  strain  7  (both  unitless) 
are  related  by  the  constitutive  relations 
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Ari  1 

oq N  — b  Oi2q~ 
t  t 


(4.4) 
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Figure  4.2  State  of  Stress  of  Differential  Wall  Element  [17] 


7  -  otiN—  +  (*4 q- 
£  £ 


(4.5) 


where  t  is  the  wall  thickness,  and  ot\,  a2,  and  0:4  are  elastic  constants  yet  to  be  defined. 
We  can  solve  for  N  and  7  by  rearranging  Eqns.  (4.4)  and  (4.5)  to  get 


N  =  fae  t  +  foq 


(4.6) 


7  =  +  Pi  q- 


(4.7) 


where  Pi  are  given  by 


(4.8) 
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(4.9) 


y04  =  «4  -  —  (4-10) 

ai 

To  determine  the  elastic  constants  07,  0:2,  and  04  we  need  to  write  the  full  anisotropic 
strain-stress  (constitutive)  equations  for  plane  stress  (<733  =  0): 

e  =  SuN±  +  S12N,j  +  Suqj  (4.11) 

e,  =  S12Nj  +  S22N±  +  Sml-  (4.12) 

7  =  SUN^  +  S2iNsl-  +  Suq\  (4.13) 

where  S  is  the  matrix  inverse  of  the  elasticity  matrix,  called  the  compliance  matrix.  By 
ignoring  the  transverse  strain  es,  setting  the  transverse  tension  flow  Ns  to  zero,  and  com¬ 
paring  with  Eqns.  (4.4)  and  (4.5)  we  get 


« 1  =  Sn 


(4-14) 


«2  =  S14 


(4.15) 
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Q?4  =  S44 


(4.16) 


Therefore,  by  determining  the  compliance  matrix  S  for  the  given  material  we  also  determine 
the  coupling  elastic  constants  a 2,  and  /?2  • 

4.  1.3  Analysis. 

4. 1.3.1  Preliminary  Considerations.  We  can  establish  the  differential 
equations  of  equilibrium  for  the  differential  element 

p-  =  0  (4.17) 

az 


dq 

ds 


dN  n 

+  a7“° 


(4.18) 


From  Eqn.  (4.17)  we  can  see  that  the  shear  flow  q(s)  is  a  function  of  s  only.  Integrating 
Eqn.  (4.18)  we  get 


rs  dN 

Q(s)=qo-J ^  ~Q^dx  (4-19) 

Figure  (4.3)  represents  the  differential  cross  section  dz  of  the  single  cell  beam,  indicating 
the  longitudinal  shear  flow  qo  at  s  =  0,  and  the  differential  change  in  the  tension  flow 
N  with  dz .  In  the  classical  theory  of  isotropic,  homogeneous,  thin-walled  beams,  the 
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cross-sectional  normal  stress  is  assumed  to  vary  linearly  with  x  and  y  [17].  That  is, 

e  =  e0  -  yKx  -  xKy  (4.20) 

where  eo,  Kx,  and  Ky  are  functions  of  z.  This  formulation  assumes  only  linear  midplane 
translational  strain  eo  in  the  2  coordinate  direction,  as  well  as  linear  curvatures  (or  also 
called  rotations  in  some  literature)  Kx  and  Ky  of  the  cross  section  at  z  around  the  x  and  y 
coordinate  axes,  respectively.  In  this  way  the  warping  functions  are  ignored  from  the  dis¬ 
placement  formulation  that  will  lead  to  less  accurate  answer  for  the  torsional  displacement. 
See  reference  [1]  for  accounting  for  bi-quadratic  warping  functions  and  the  qualitative  im¬ 
provement  they  yield  in  the  results.  Let’s  substitute  Eqn.  (4.20)  into  Eqn.  (4.6)  and 
differentiate  with  respect  to  z  to  get 

N  =  fat{eo  -  yKx  -  xKy)  +  foq  (4-21) 
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—  =  /31t(e,0-yK'x-xK'y)  (4.22) 

Substituting  Eqn.  (4.22)  into  Eqn.  (4.19)  we  get 

q(s)  =  q0  -  [  Pit(e'0  -  yK'x  -  xK'y)ds  (4.23) 

Jo 

We  can  simplify  Eqn.  (4.23)  by  declaring 

q(s)  =  qo-  eoai(s)  +  K'xa2{s )  +  K'ya3{s)  (4.24) 


where 


L  W  =  / 
Jo 


ai(s)  =  /  Pitds 


l(s)=  [ 
Jo 


a2{s)  =  /  yfiitds 


03(5)  =  [ 
Jo 


xf$\t  ds 


(4.25) 


(4.26) 


(4.27) 


Equations  4.25  through  4.27  are  solved  individually  and  consecutively  for  all  four 
surfaces  of  the  trapezoid,  yielding  four  equations  for  each  a*(s)  equation.  Each  surface 
will  have  its  own  contribution  that  is  added  to  the  contribution  of  all  surfaces  previously 
evaluated.  The  method  of  solving  Eqns.  (4.25)  through  (4.27)  is  demonstrated  via  ai(s). 
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The  solution  of  a2(s)  and  03(5)  follows  the  same  technique. 


ai(S)4=  r 
Jo 


Pitsds 


(4.28) 


r2dT  rs 

ai(s)r  =  /  Pitsds+  /  fifads 

Jo  J2d 


.  (4.29) 


n&T 


r2dT  r2dT-\-h2  rs 

ai(s)b  =  /  f3itsds  +  /  fiihds  +  /  Pi  tsds  (4.30) 

Jo  JldT  J2d 


l2dT+h2 


r2  dT  r2dT+h2  r4dT~{-h2  rs 

ai(s)^  =  /  f3itsds  +  /  /Ji^ds  +  /  /?i£sds  +  /  /3i£ids  (4.31) 
Jo  J2dT  J2dT+h2  J  4dT+h2 

where  £s,  ti,  and  £2  are  the  thickness  of  the  top  and  bottom  skin,  the  main  spar  (left  side) 
and  rear  spar  (right  side),  respectively.  The  superscripts  r,  6,  and  l  of  a^(s)  are  for 
the  top,  right,  bottom,  and  left  sides  respectively.  To  solve  Eqns.  (4.28)  through  (4.31)  I 
suggest  the  use  of  any  symbolic  solver  algorithms  currently  available,  such  as  MathCad, 
Matlab  or  Mathematica.  The  symbolic  results  of  Eqns.  (4.28)  through  (4.31)  are  relatively 
simple  and  will  not  be  reproduced  here. 

1.3.2  Equations  of  Static  Equivalence .  At  any  cross  section  of  the 
beam,  the  longitudinal  force  per  unit  length  N  must  be  statically  equivalent  to  the  cross- 
sectional  extension  force  P,  the  cross  sectional  bending  moment  about  the  x  axis  Mx,  and 
the  cross  sectional  bending  moment  about  the  y  axis  My  [17].  By  using  Eqn.  (4.21)  to 
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eliminate  N  from  the  integrals,  we  get 


P  =  <j>  Nds  =  j>  -  yKx  -  xKy)  +  ^2q]ds 


(4.32) 


Mo 


=  -  j)  yNds  =  -  j>  y\Pit(eo  -  yKx  -  xKy )  +  hq}ds  (4.33) 


My  =  -  j>  xNds  =  -  j)  x[Pit(eo  -  yKx  -  xKy)  +  hq)d& 


(4.34) 


In  order  to  express  these  equations  in  a  matrix  form,  let’s  first  write  out  Eqns.  (4.32) 
through  (4.34): 


eo  Pitds  -  Kx  <j>  yfiitds  -  Ky  j)  xfiitds 
-eo  j)  yPitds  +  Kx  j)  y2/3itds  +  Ky  j)  xy(3\tds 
-eo  xfiitds  +  Kx  j)  xyf3\tds  +  Ky  j)  x20itds 


P  -  foqds 
Mx  +  yfoqds 


=  My  +  f  xfcqds 


/■ 


(4.35) 

(4.36) 

(4.37) 


After  multiplying  the  last  two  equations  by  —1,  we  can  write  them  in  a  more  manageable 
matrix  form 


hi  h2  hz 

f  \ 

eo 

r  \ 

P-Qi 

hi  hi  hz 

< 

Kx 

►  =  < 

1 

1 

<o 

to 

hi  h2  hz 

Ky 
\  > 

^  ~My  ~  Qs  J 

(4.38) 
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where  the  b  matrix  is  made  up  of  the  following  elements: 


hi  =  §  Pitds  bi2  =  ~  §  yP lids  613  =  -  §  xfii tds 

hi  =  j>  yftitds  b22  -  -  §  y2Pitds  623  =  -  §  xyfiitds  (4-39) 

631  =  §  xp\ tds  632  =  -  f  xyPitds  633  =  -  f  x2/3itds 

We  can  certainly  solve  [btj]  by  hand,  but  they  are  best  left  to  computer  symbolic  solvers. 

It  is  important  that  we  observe  these  integrals  are  full  path  integrals,  and  the  limits  of 
integrations  are  as  defined  in  Eqn.  (4.3).  Due  to  the  asymmetry  about  the  y  axis  (we 
defined  the  cross  section  of  the  torque  box  as  a  trapezoid) ,  the  b  matrix  is  not  diagonal. 
The  off-diagonal  elements  are  all  zero,  except  613  =  —631  ^  0.  Also,  to  complete  the 
solution  of  Eqn.  (4.38),  we  have 

Qi  =  §  qfods  ,  Q2  =  f  yfoqds  ,  Qz  =  §  xfaqds  (4-40) 

However,  for  the  homogeneous,  isotropic  case,  Qi  =  0,  because  the  coupling  elastic  constant 
«2  =  0;  therefore,  fh  =  0.  In  the  non- homogeneous,  anisotropic  case  discussed  later  (see 
Section  4.2),  the  S  compliance  matrix  is  fully  populated,  and  a2  /  0;  therefore,  7^  0. 

In  order  to  solve  for  the  unknown  functions  eo  (longitudinal  strain  in  the  2  direction), 
Kx  (rotation  about  the  x  axis)  ,  and  Ky  (rotation  about  the  y  axis)  of  Eqn.  (4.20)  we  invert 
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[bij]  to  get  [ciij].  That  is 


r  y 

- 

- 

f  \ 

eo 

on 

012 

013 

P-Qi 

Kx 

i  = 

021 

022 

«23 

< 

—Mx  —  Q2  > 

Ky 

v  y ) 

«31 

032 

“33 

~My  ~  Q3 

< 

(4.41) 


Differentiating  Eqn.  (4.41)  with  respect  to  z  (P  —  const ,  dMx/dz  =  —Vx  =  const, 


dMy/dz  =  —  Vy  =  const),  we  get 


\ 

“ 

“ 

( 

eo 

Oil 

012 

013 

0 

K 

>  = 

021 

022 

023 

i 

> 

Ky 

V  j 

031 

032 

033 

Vx 

W  J 

(4.42) 


The  shear  flows  q(s)  must  be  statically  equivalent  to  the  applied  moment  M  [20]. 
Here  we  define  the  quantity  p(s)  as  the  perpendicular  distance  from  the  origin  to  the 
point  on  the  cross  section  defined  by  the  path  coordinate  s.  In  case  of  the  trapezoid  cross 
section,  p(s)  =  const  for  all  sides.  For  the  left  and  right  sides  p(s)  =  d.  For  the  top  and 
bottom  sides,  p(s)  can  be  very  closely  approximated  by  p(s)  =  |(hi  +  ^2)-  Because  p(s) 
is  a  measure  of  distance,  p(s)  =  d  for  the  left  side  as  well  as  for  the  right  side.  Using 
Eqn.  (4.24)  the  expression  for  the  applied  moment  can  is  written  as 


M  =  <j>  p{s)  [90  -  e'oOi(s)  +  K'xa2(s)  +  ^o3(s)]  ds 


(4.43) 
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from  which  we  can  directly  solve  for  qo  by 


p{s)qQds  =  M  +  ef0 


p{s)ai{s)ds )  -K'x 


p(s)a2(s)dsSj  -  K'y  (J>  p(s)a3(s)ds 

(4.44) 


*7°  ^p(s)ds  ^  6°a4  KXa5  Kyaf>\ 


(4.45) 


where  the  constants  as,  and  a§  are  given  by 


a4  =  j>  p(s)ai(s)d6 


(4.46) 


a  5 


=  j)  p(s)a2(s)di 


(4.47) 


j>  p(s)a3(s)ds 


(4.48) 


The  quantity  j^yf$  in  Eqn.  (4.45)  is  the  reciprocal  of  twice  the  torque  box  area,  that  is 


^p(s)ds  2A 


(4.49) 


4. 1.4  Rate  of  Twist.  The  classical  theory  of  thin- walled  beams  relate  the 

rate  of  twist  of  a  cell  to  the  shear  strains  in  the  wall  of  the  cell  [26].  Because  we  are 
determining  the  static  deformation  of  the  section,  the  rate  refers  to  the  change  of  the  angle 
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along  the  longitudinal  coordinate  as  opposed  to  the  change  of  angle  in  time.  We  write 
the  relation  as 


d^ 
d  z 


(4.50) 


where  <f>(s)  is  the  angle  of  twist.  In  order  to  define  the  shear  strain  7(5),  we  need  to  revisit 
Eqn.  (4.7),  reproduced  here  for  convenience: 


7  =  -fo  +  (4-51) 

The  longitudinal  strain  e  was  defined  by  Eqn.  (4.20),  and  the  shear  flow  q(s)  was  given  in 
Eqn.  (4.24).  Using  these  two  equations  in  Eqn.  (4.7),  we  have 

7(s)  =  ~/32{e o  -  yKx  -  xKy)  +  /34  [q0  -  eoai(s)  +  K'xa2{s )  +  K'yaz(s)]  1  (4.52) 

Substituting  Eqn.  (4.52)  into  Eqn.  (4.50),  we  get 

^  [—e0cn  +  Kxc 2i  +  Kyc 3i  +  q^di  +  e’0d2  +  K^d^  +  Kydi]  (4.53) 

where  the  constants  <i,  are  given  by 


di 


(4.54) 


d2  —  — 


(4.55) 
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d3  =  a2(s)P4-ds 


(4.56) 


C?4 


a3(s)/34-ds 


(4.57) 


The  methodology  solving  Eqns.  (4.54)  through  (4.57)  follows  the  technique  outlined  in 
Eqns.  (4.28)  through  (4.31).  For  example,  solving  for  d2  for  each  side  of  the  trapezoid,  we 
must  use  the  corresponding  expression  for  ai(s)  from  Eqns.  (4.28)  through  (4.31),  and  the 
corresponding  value  for  the  wall  thickness  t(s).  That  is 


d2  = 


p2dT  i 

/  a[(s)/34-ds 

Jo  ZS 

r2dT+h2  i 

/  a^&j-ds 

J2dT 


h 


r4dT-\-h2  i 

/  a?(s)Ai-ds 

J2dT+h2  Zs 


L 


4dT+/i2+/ii  | 

a[(s)Pi—ds 
1 1 


4dT+/i2 


(4.58) 


The  constants  Cij  in  Eqn.  (4.53)  are  all  identically  zero  for  the  isotropic  case,  because  the 
coupling  elastic  constant  /3 2  is  zero.  Therefore,  the  equation  for  the  rate  of  twist  for  the 
single-cell  isotropic  thin-walled  beam  can  be  simplified  from  Eqn.  (4.53)  as 


=  ttt  Mi  +  e'o d2  +  K'xdz  +  K'ydA) 


(4.59) 


where  the  constants  di  are  known  through  Eqns.  (4.54)  through  (4.57),  the  constant  shear 
flow  go  is  given  by  Eqn.  (4.45),  and  the  constants  e'0,  K'x,  and  K'y  are  found  from  Eqn.  (4.42). 
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4- 1-5  Isotropic  Case  Identified.  The  motivation  behind  developing  the 
homogeneous,  isotropic  solution  using  the  full  anisotropic  methodology  was  to  be  able 
to  check  the  result  and  accuracy  of  the  solution.  The  answer  from  the  single-cell  tor¬ 
sion  solution  (see  Chapter  III)  will  serve  as  a  ’yardstick’  for  the  result  obtained  from  the 
homogeneous,  single-cell  torsion  solution. 


4.2  N on-homogeneous  Anisotropic  Single  Cell  Beam 


4*2.1  Considerations.  The  definition  of  non-homogeneous  cross-section  in 
this  context  means  that  the  coupling  elastic  constant  /?2 ,  and  therefore  a2  are  functions 
of  the  path  length  coordinate  s.  Due  to  the  non-zero  (and  also  not  90  degrees)  composite 
fiber  orientation  angle  0,  the  elastic  coupling  constant  will  change  sign  as  the  fiber 
angle  changes  sign  with  respect  to  the  z  structural  axis.  That  is,  the  angle  6  between  the 
structural  axis  z  and  the  principal  fiber  direction  is  greater  than  zero  if  the  fiber  is  rotated 
with  respect  to  the  z  axis  in  a  positive  sense  in  the  right-handed  coordinate  system.  The 
positive  rotation  is  defined  as  the  direction  of  the  vector  cross  product  of  these  two  axes 
that  define  the  plane  of  the  rotation  angle  6.  If  the  fiber  is  wound  continuously  around  the 
box  beam  we  have: 


>0  : 

So  <  S  <  Sl 

>0  : 

Si  <  s  <  S2 

<0  : 

S2  <  s  <  s 3 

<0  : 

S3  <  S  <  S4 

(4.60) 


where  the  limits  Si  are  identified  by  Eqn.  (4.3). 
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4-2.2  Analysis.  Because  the  method  of  analysis  outlined  in  Section  4.1  is 
derived  for  a  fully  anisotropic  case,  then  applied  to  an  isotropic  example,  all  the  steps 
outlined  in  that  section  apply  to  the  current,  non-homogeneous,  anisotropic  solution.  In 
Eqn.  (4.40)  we  found  the  integrals  that  define  Qx,  Q2,  and  Q3.  In  this  case,  these  constants 
are  not  zero,  because  the  elastic  coupling  constant  /?2  is  not  equal  to  zero.  This  comes  from 
the  fact  that  the  compliance  matrix  S  is  fully  populated;  therefore,  a2  =  Si 4  7^  0.  Because 
a\  is  always  non-zero,  using  Eqn.  (4.9),  the  constant  /%  is  defined  and  non-zero. 

Recall  that  up  until  Eqn.  (4.52),  none  of  the  formulas  involve  /?2j  therefore,  the  results 
obtained  in  Section  4.1  remain  valid.  From  Eqn.  (4.40)  we  have 


qo§ P2ds-e'Q§  ax{s)P2ds  + K'x§  a2{s)P2ds  + K'y§  a3(s)/52ds  =  Qi 

Qo  $ yfads  -e'0§ yai{s)fi2ds  +  K'X§  ya2(s)/32ds  +  K'y§ ya3(s)/32ds  =  Q2 

q0  §  xf32ds  —  e'0  §  xai(s)f32ds  +  K'x  § xa2(s)/32ds  +  K'y  $ xa,3(s)/32ds  =  Qz 


(4.61) 


which  we  can  write  in  matrix  notation  as 


f  ^ 

- 

- 

go 

f  \ 

Cll 

C12 

C13 

C14 

e0 

Qi 

C21 

C22 

C23 

C24 

< 

K 

>  =  i 

Q2  1 

C31 

C32 

C33 

C34 

Qz 

K'y 
\  U  J 

V 

(4.62) 
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where  the  elements  of  the  c  matrix  are  now  defined  as 


cn  =  f  /32ds  C12  =  -  §  ai(s)/32ds  C13  =  §  a2(s)/32ds  cu  =  /«3(s)^2ds 

c2i  =  § y/32ds  C22  —  ~  § yai(s)p2&s  c2 3  =  § ya2(s)/32ds  C24  =  / ya3(s)/%ds  (4.63) 

C31  =  j>  xfyds  C32  =  -  f  xai(s)/02ds  c33  =  <f  xa2(s)/3 2ds  C34  =  § xa3(s)(32ds 

Finding  Cjj  by  the  path  integrals  requires  some  clarification.  For  example,  C23  is  summed 
clockwise  around  the  perimeter,  starting  from  s  —  0.  For  each  consecutive  side  of  the 
trapezoid,  we  must  use  the  corresponding  expression  for  a2(s).  Therefore,  the  constant  c23 
is  given  by 


C23 


-L 


31  f/l2  hl(s-Td)  +  \(h2  +  h  4) 


4  dT 


^2a2(s)ds  + 


+  f  (2Td  -  s  +  ^/i2)/32G2(s)ds  + 

V  Si  w 

+  /  _  4dT  1  _  s  +  ^2)  +  4(^2  +  hi) 

+  f  (s-  4Td  -h2-  ^/ii)(-/?2)4(s)ds 


(-/02)o2(s)ds  + 


(4.64) 


where  the  limits  of  integrations  s*  are  given  by  Eqn.  (4.3),  and  the  superscripts  on  d2(s) 
refer  to  the  top,  right,  bottom,  and  left  sides  of  the  trapezoid,  respectively.  Observe  that 
the  sign  of  /?2  changes  on  the  bottom  and  the  left  sides  of  the  trapezoid.  This  is  the  effect 
of  the  non- homogeneous  construction,  when  the  bottom  side  is  the  mirror  image  of  the 
top  side,  that  is,  the  fiber  orientation  angle  on  the  bottom  surface  (with  respect  to  the 
right-handed  structural  coordinate  system)  is  equal  and  opposite  to  that  of  the  top  surface 
[17].  For  the  same  reason,  the  sign  of  /?2  changes  from  the  right  side  to  the  left  side.  The 
solution  for  the  rest  of  the  c  matrix  follows  the  same  procedure,  without  exceptions. 
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4-2.3  Rate  of  Twist. 


The  rate  of  twist  of  the  homogeneous  single-cell 


section  was  developed  in  Section  4.1.4.  Equations  (4.50)  and  (4.53)  remain  valid  in  this 
case,  except  now  we  must  account  for  the  constants  Cy  as  well.  Because  the  functions  Qt 
are  not  zero,  we  first  need  to  solve  for  eo,  Kx.  and  Ky.  From  Eqn.  (4.41)  we  have 


" 

“ 

'I 

eo 

Oil 

012 

013 

1 

<0 

Kx 

>  = 

021 

022 

023 

< 

-Mx  -  Q2  \ 

Ky 

031 

032 

033 

s - 

l 

1 

CO 

(4.65) 


where  Qi  are  given  by  Eqn.  (4.62).  Once  we  perform  the  matrix  multiplications,  we  can 
substitute  Qi  into  Eqn.  (4.41)  and  obtain  the  expressions  for  eo,  Kx,  and  Ky.  The  final 
format  of  Eqn.  (4.65)  will  be 


eo  =  a-i\{P  ~  Qi)  ~  ai2 {Mx  +  Q2)  —  ais(My  +  Q 3)  (4.66) 


Kx  —  a2i(P  —  Qi)  —  a  +  Q2 )  —  a2s{My  +  Q  3)  (4-67) 

Kx  =  031  {P  ~  Qi)  —  032 (Mx  +  Q2)  —  asz(My  +  Q3)  (4.68) 

Now  that  eo,  KX)  and  Ky  are  given  by  Eqns.  (4.66)  through  (4.68),  e'0,  K'x,  and  K'y  are 
determined  by  Eqn.  (4.42),  the  constants  Cy  are  given  by  Eqn.  (4.63),  the  constants  di 
are  known  by  Eqns.  (4.54)  through  (4.57),  we  can  write  the  rate  of  twist  for  the  non- 
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homogeneous,  anisotropic,  single-cell  section  as 


d  z 


1 

2  A 


(— eocn  +  KxC2i  +  KyCn  +  <70^1  +  ef0d2  +  Kfxds  +  Kfyd 4) 


(4.69) 


4*3  Discussion 

To  evaluate  Eqns.  (4.8)  through  (4.69)  simultaneously  for  the  given  geometry  is 
a  formidable  as  well  as  lengthy  task.  The  evaluation  of  the  full  and  the  partial  path 
integrals  could  easier  be  done  by  using  symbolic  solvers.  The  symbolic  results  of  Eqns.  (4.8) 
through  (4.69)  were  entered  into  a  Matlab  subroutine  called  ANISOTORSION.COMPZ 
that  calculates  the  tip  twist  angle  of  the  non-homogeneous  (see  Section  4.2)  anisotropic, 
single-cell,  closed-section  box  beam.  This  code  can  also  be  used  to  calculate  the  tip  twist 
of  a  homogeneous  anisotropic  beam,  since  in  this  case  the  coupling  constant  /?2  is  zero  (see 
Section  4. 1.3. 2). 

Another  Matlab  subroutine  called  LAYUP  [23]  was  applied  in  order  to  calculate 
the  engineering  properties  of  the  fully  anisotropic  composite-laminate,  single-cell,  closed 
section  box  beam,  so  that  the  stiffness  and  compliance  matrices  of  the  laminae  in  the  1-2 
principal  axes  can  be  transformed  into  the  x-y-z  structural  axes  of  the  laminate,  and  sent 
to  the  subroutine  ANISOTORSION.COMPZ. 

A  driver  program  called  COMPL_ANISO_COPMZ  was  written  in  order  to  define  the 
material  properties  in  the  1-2  principal  axes  of  the  anisotropic  composite  lamina  materials 
of  choice.  All  three  subroutines  and  programs  are  listed  in  Appendix  C. 
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V .  Induced  Strain  Actuation 


5.1  Single- Cell  Beam  with  PZT  Lamina 

5.1.1  PZT- Composite  Layup.  In  order  to  produce  induced  strain  in  the 
top  and  the  bottom  surfaces  of  the  single-cell  box  beam,  a  single  layer  of  continuous-sheet 
PZT  lamina  is  incorporated  in  the  composite  laminate  of  the  top  and  the  bottom  surfaces 
of  the  beam  as  shown  in  Figure  (5.1).  In  addition  to  the  assumptions  made  in  Section  4.1.1, 
it  is  assumed  that  the  PZT  layer  is  perfectly  bonded  within  the  host  structure,  that  is  no 
slip,  disbond,  or  shear-lag  is  accounted  for.  The  PZT  lamina,  when  used  as  an  actuator 


Figure  5.1  Torque  Box  with  Embedded  Piezoelectric  Lamina 

under  tension,  pulls  the  host  structure  inward;  when  under  compression,  it  pushes  the 
host  structure  outward.  If  the  top  PZT  lamina  is  subjected  to  an  electric  field  that  creates 
compression,  while  the  bottom  PZT  lamina  experiences  tension,  the  entire  single-cell  beam 
will  be  subjected  to  a  pitching  moment  about  the  negative  z  structural  axis.  This  way 
the  beam  is  torqued  nose  up,  and  the  angle-of-attack  increases.  If  the  electric  fields  are 
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switched  on  the  PZT  laminae,  the  moment  changes  direction,  and  the  beam  is  subjected 
to  a  positive  torque  that  pitches  the  beam  nose  down.  See  Figure  (5.2)  for  a  simplified 
representation  of  the  PZT  actuation  of  a  generic  (geometry  is  not  related  to  the  simplified 
trapezoid  section)  host  structure. 

Keajelectric 

Actuator 


Figure  5.2  Piezoelectric  Lamina  Configuration  [12] 


To  achieve  high  PZT  actuator  effectiveness,  a  PZT  lamina  with  high  piezoelectric 
constant  (strain  coefficient)  d 33  must  be  selected.  If  d^s  of  the  PZT  is  small,  a  large 
voltage  is  required  to  produce  strain  in  the  PZT.  If  c/33  is  large,  a  small  amount  of  voltage 
is  sufficient  to  produce  the  necessary  strain  [9].  The  PZT  must  also  have  high  Young’s 
modulus  of  elasticity  E  compared  to  the  host  structure,  so  that  a  large  fraction  of  the 
strain  produced  by  the  electric  field  can  be  transferred  to  the  host  structure  [9]. 

Piezo-fiber  composites  are  becoming  more  readily  available,  as  well  as  less  expensive. 
For  this  study,  two  PZT  fibers  were  considered,  with  almost  identical  material  properties. 
The  first  PZT  lamina  selected  is  an  Active  Fiber  Composite  (AFC)  PZT  5H  fibers  with 
thermosetting  epoxy  resin  matrix  and  etched  copper/Kapton  interdigitated  electrodes.  A 
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Poling  Direction 


Figure  5.3  Diagram  and  Cross-section  of  the  AFC  Lamina  [5] 


Table  5.1  PZT-Fiber  Composites 


Fiber  Type 

Ei 

(Msi) 

e2 

(Msi) 

Gu 

(Msi) 

^12 

t 

(in) 

^33 

(pm/V) 

dn 

(pm/V) 

AFC  Lamina 

4.6786 

2.4173 

0.5802 

0.30 

0.0065 

180 

-50 

G-1195  PZT 

5.4389 

2.0305 

0.5511 

0.30 

0.0080 

200  -  400 

0.1^33 

representative  example  of  the  AFC  PZT  is  illustrated  in  Figure  (5.3).  The  other  is  a 
generic  piezo-fiber  composite  from  reference  [3].  The  material  properties  of  the  two  PZT 
composites  are  summarized  in  Table  (5.1). 

Before  attempting  to  derive  the  torsional  formula  for  this  PZT-composite  beam, 
we  have  to  incorporate  the  PZT  engineering  properties  into  the  those  of  the  composite 
laminate  alone.  We  achieve  this  by  modifying  the  Matlab  routine  used  for  calculating  the 
engineering  properties  of  the  laminate  alone  so  that  it  accounts  for  the  single  ply  of  PZT. 
For  sake  of  simplicity,  but  without  any  loss  of  generality,  the  PZT  lamina  was  applied  to 
the  center  of  the  laminate  by  replacing  the  middle  composite  layer. 


5.1.2  Constitutive  Relations.  The  plane  stress  constitutive  relations 
(strain-stress)  of  Section  4.1.2  have  to  be  modified  to  account  for  the  PZT  layer,  that 
may  not  have  the  same  material  properties  as  the  composite  lamina.  Therefore,  we  write 
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Eqns.  (4.4)  through  (4.5)  with  the  PZT  strain  (ep),  and  the  PZT  shear  (7^)  in  the  material 
principle  directions 


e  —  Oi\N-  +  Cp 

u  £ 


(5.1) 


7  =  a2N-  +  a4<?-  +  7P 

Using  the  same  argument  preceeding  Eqns.  (4.6)  and  (4.7)  we  have 


N  =  P\e  t  +  fa q  -  Piep  t 


(5.2) 


(5.3) 


7  —  —  fat  +  +  @2ep  +  Ip  (5-4) 

where  all  the  constants  evaluated  to  be  the  same  as  defined  in  Eqns.  (4.8)  through  (4.10). 
Writing  the  full  anisotropic  strain-stress  (constitutive)  equations  for  plane  stress  and  com¬ 
paring  them  to  Eqn.  (5.1),  we  obtain  the  same  results  as  in  Eqns.  (4.14)  and  (4.16). 
Therefore,  the  addition  of  the  PZT  term  in  the  Eqns.  (5.1)  and  (5.2)  did  not  change  our 
coupling  coefficients  defined  previously  in  Eqns.  (4.8)  and  (4.10). 

5.1.3  Analysis. 

5. 1.3.1  Preliminary  Considerations.  By  establishing  the  same  differ¬ 
ential  equations  of  equilibrium  considered  in  Section  4. 1.3.1,  and  using  Eqns.  (4.17)  and 
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(4.20),  we  obtain  the  expression  for  the  force  per  unit  length  N,  namely 


N  =  /3it(e0  -  yKx  -  xKy)  +  fcq  -  P\ep  t 


(5.5) 


Differentiating  with  respect  to  z,  we  get 


dN 

dz 


=  /31t(e'0-yK'x-xK'y)-l3le'pt 


(5.6) 


where  the  derivative  of  the  piezoelectric  strain  with  respect  to  the  spanwise  coordinate  z 
is  zero  as  long  as  uniform  voltage  is  applied  to  the  PZT  lamina.  That  is 

4-&-0  M 

With  this  assumption,  we  obtain  the  expression  for  dN/dz 

TP  =  PM  -  yK'x  -  xK'y)  (5.8) 


that  is  identical  to  the  result  derived  in  Section  4. 1.3.1,  Eqn.  (4.22).  Therefore,  the  rest  of 
the  derivation  for  Eqns.  (4.24)  through  (4.31)  is  as  outlined  before. 

5. 1.3. 2  Equations  of  Static  Equivalence.  The  cross  sectional  loadings 
will  now  have  to  be  modified  to  account  for  strain  actuation  as  well.  Equations  (4.32) 
through  (4.34)  are  written  in  a  form 


/ 


[/M(e o  -  yKx 


-  xKy)  +  /?2 q  ~  j3iept]ds 


(5.9) 
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(5.10) 


Mx  =  -  j)  yNds  -  -  j>  y[Pit(e0  -  yKx  -  xKy)  +  /3 2q  ~  /3iept\ds 


M, 


=  -  j)  xNds  =  -  j)  x[/3it(eo  -  yKx  —  xKy)  +  fa Q  —  faept]ds  (5.11) 


that  can  be  written  as  (see  also  Eqn.  (4.38)) 


1 

00 

*0 

CM 

-O 

rH 

1 - 

/  \ 

€0 

f  \ 

P  —  Qi  +  Z\ 

hi  hi  ^23 

Kx 

i  =  < 

—Mx  —  Q2  +  Z2  > 

hi  h2  633 

Ky 

\  j  ) 

—  My  —  Q%  +  Z3 

k  / 

(5.12) 


where  the  [5]  matrix  is  given  by  Eqn.  (4.39),  and  the  quantities  Z{  are  given  by 


Z\  =  §  faeptds  Z2  =  f  yfaeptds  Z$  =  <f  xfaeptds  (5.13) 

The  full  path  integrals  of  Eqn.  (5.13)  must  be  evaluated  the  same  way  we  treated  the  path 
integral  in  Eqn.  (4.58).  In  doing  so  we  find 


Zi  =  4fatsepTd  Z2  =  0  Z3  =  0  (5-14) 

Evaluating  the  derivatives  of  Eqn.  (5.14)  with  respect  to  z  we  arrive  to 


z[  =  0  z'2  =  d  z’3  =  0 


(5.15) 
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After  rearranging  Eqn.  (5.12),  and  differentiating,  we  get 


f  > 

- 

“ 

/  > 

e0 

an 

012 

Ol3 

0 

< 

>  = 

021 

022 

023 

< 

Vy  > 

031 

O32 

033 

14 

^  / 

(5.16) 


which  is  identical  to  Eqn.  (4.42). 

5.1.4  Rate  of  Twist.  In  light  of  the  results  of  Section  5.1.3,  all  equations 
apply  until  Eqn.  (4.50),  which  is  reproduced  here  as 


=  U  /7(s)di 


Cty  1 
dz 


(5.17) 


where  the  shear  strain  is  given  by 


7  =  a2A^  +  a^qj  +  7p 


(5.18) 


This  can  also  be  written  as 


7  =  -fat  +  faqj  +  fa  ep  +  7 p 


(5.19) 


which  is  Eqn.  (5.4)  from  Section  5.1.2.  Substituting  from  Eqns.  (4.20)  and  (4.24),  we  get 


700  =  ~fa(eo  ~  yKx  -  xKy )  +  /?4  [q0  -  e'0ai(s)  +  AT'a2(s)  +  K'ya3{s )]  -  +  faep  +  7p 

(5.20) 
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which  modifies  Eqn.  (4.53)  to 


d0 
d  z 


1 

2 A 


-eocn  +  KxC2i  +  KyC^i  +  qodi  +  e'0d2  +  K'xdz  +  K'ydA  +  c\\  j)  epds  +  j)  7pd s 


where  the  constants  d{  are  as  given  by  Eqns.  (4.54)  through  (4.57).  The  full  path  integral  of 
the  piezoelectric  strain  yp  is  evaluated  the  way  it  was  demonstrated  for  C23  of  Eqn.  (4.64). 
Assuming  PZT  lamina  only  in  the  skins,  we  obtain 


7p(s)ds  =  47 pTd 
j)  ep(s)ds  = 


4  evTd 


(5.22) 

(5.23) 


where  ep  and  7 p  is  the  PZT  strain  and  shear  in  the  structural  axes,  respectively.  For  a 
homogeneous  cross  section,  all  Cij  are  zero,  therefore 

(5.24) 


~r~  =  t"T  Ml  +  t'0d2  +  K'xd3  +  K'dA  +  4 lpTd) 


5.2  N on-homogeneous  PZT- Composite  Beam 

5.2.1  Analysis.  For  a  non- homogeneous  cross  section,  we  follow  Section  4.2. 
Because  all  Qi  s  are  still  defined  in  the  PZT  composite  case  the  same  way  they  were  in  the 
composite-only  case,  Eqns.  (4.60)  through  (4.63)  still  hold. 

5.2.2  Rate  of  Twist.  Recall  that  for  the  non- homogeneous  cross  section, 
^  0.  However,  Eqn.  (4.65)  is  now  extended  with  the  summation  of  the  longitudinal 
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strain  terms,  Z%.  That  is, 


Oil 

«12 

021 

022 

031 

032 

P-Qi  +  Zi 
—Mx  —  Q2  +  Z2 

—  My  —  Qz  +  Z3 


(5.25) 


where  Z%  are  given  by  Eqn.  (5.13).  Using  the  results  of  Eqn.  (5.14)  we  can  evaluate 
Eqn.  (5.25)  and  get 


eo  =  ai  i(P  —  Qi  +  %i)  —  &12  (Mr  +  Q2)  —  ^i3  (My  +  Qz)  (5.26) 


Kx  =  Q2i{P  —  Qi  +  Z\)  ~  022 {Mx  +  Q2)  —  023 {My  +  Qz)  (5-27) 


Ky  =  azi{P  —  Qi  +  Z\)  —  az2(Mx  +  Q2 )  —  azz{My  +  Qz )  (5.28) 

where  Z2  =  Z3  =  0  as  per  Eqn.  (5.14).  Now  that  eo,  Kx,  and  Ky  are  given  by  Eqns.  (5.26) 
through  (5.28),  e'0,  K'x,  and  K'y  are  determined  by  Eqn.  (5.16),  the  constants  are  given 
by  Eqn.  (4.63),  the  constants  di  are  known  by  Eqns.  (4.54)  through  (4.57),  we  can  write 
the  rate  of  twist  for  the  non- homogeneous,  anisotropic,  single-cell  section  as 


(5.29) 
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5.3  Determining  the  PZT  Strain  Tensor 


Using  Eqn.  (5.29)  we  can  determine  the  angle  of  twist  due  to  PZT  actuation  of  the 
non- homogeneous,  anisotropic,  single-cell  beam  with  PZT  actuator  lamina  embedded  in 
the  laminate  center  of  the  top  and  bottom  skins.  All  variables  and  coefficients  of  Eqn.  (5.29) 
are  known  through  the  derivation  presented  in  Chapters  IV  and  V,  except  the  PZT  strain 
ep,  and  the  PZT  shear  7p.  Before  calculating  the  twist  angle  due  to  an  applied  electric 
field  (Voltage),  we  need  to  find  ep  and  jp  of  the  PZT.  Because  the  PZT  actuator  lamina  is 
at  an  angle  with  respect  to  the  material  axes  of  the  torquebox,  we  first  need  to  determine 
ep  and  7P  in  the  PZT  lamina  principal  axes,  then  transform  these  strain  and  shear  to  the 
structural  axes. 

For  the  sake  of  preserving  the  conventions  for  the  directions  of  piezoelectric  actuation, 
let’s  rename  our  principal  1-2-3  (fiber-transverse-out  of  plane)  composite-fiber  material 
directions  to  3-1-2  (fiber-transverse-out  of  plane)  composite-PZT  material  directions,  as 
demonstrated  in  Figure  (5.4).  This  way  we  can  retain  the  subscripts  on  all  the  piezoelectric, 

and  electric  field  coefficients  to  be  introduced  in  the  following  section. 

Top  Etched  Copper/ 


Figure  5.4  Diagram,  and  Principal  Axes  of  the  AFC  lamina  [5] 
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5.3.1  Piezoelectric  Strain  Coefficients.  According  to  the  notation  of 
IEEE  Standard  176-1978,  the  linear,  coupled  electro-mechanical  constitutive  relations  are 


D  =  eTE  +  dT 
S  =  dtE  +  sET 


(5.30) 


where  the  independent  variable  stress  (T)  and  electric  field  (E)  and  the  dependent  vari¬ 
ables  strain  (S)  and  electric  displacement  (D)  are  related  by  the  dielectric  coefficient  ( eT ), 
piezoelectric  coefficient  (d),  and  the  compliance  matrix  (SE)  [8].  For  piezoelectric  mate¬ 
rials  of  unit  thickness  the  relation  between  actuation  strain  and  the  applied  electric  field 
(Ei,  E2,  £3  in  the  principal  1-2-3  directions  respectively)  is  given  by 


0 

0 

dzi 

0 

0 

<^31 

/  > 

Ei 

0 

0 

<^33 

e2  > 

0 

dm 

0 

£3 

diz 

0 

0 

t  > 

0 

0 

0 

(5.31) 


The  coefficients  <kn ,  c/33,  and  r/15  are  the  piezoelectric  strain  coefficients  (or  constants). 
The  constant  c/33  characterizes  the  strain  in  the  material  principal  fiber  direction,  and  the 
constant  dn  characterizes  the  strain  in  the  material  transverse  direction.  In  accordance 
with  the  assumptions  and  restrictions  made  on  the  analysis  in  Section  4.1.2,  and  to  simplify 
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the  relations  of  Eqn.  (5.31)  we  set 


E\  —  E2  —  0 


(5.32) 


Therefore,  the  PZT  strain  tensor  ei52,3  in  the  principal  directions  is  expressed  by 


el>2,3 


G?3i  0  0 

0  C?3i  0 


0  0  C?33 


(5.33) 


where  E$  is  given  by  the  applied  voltage  divided  by  the  PZT  lamina  characteristic  distance. 

If  the  piezoelectric  constants  are  of  opposite  signs,  the  applied  electric  field  creates 
extension  in  one  direction  (positive  constant),  and  contraction  (negative  constant)  in  the 
other.  If  the  constants  are  of  the  same  sign,  the  electric  field  will  induce  extension  or 
contraction  simultaneously  in  both  principal  directions.  For  the  45-degree  PZT  layup  used 
in  the  study  (see  Figure  (5.2)),  a  PZT  lamina  with  positive  0/33,  and  negative  d3 1  could 
generate  torsion  of  the  beam  more  effectively,  since  the  mechanical  strains  on  the  top 
and  bottom  surfaces  are  of  the  opposite  signs.  Which  beam  is  going  to  produce  more 
angular  twist;  however,  is  also  going  to  depend  on  the  numerical  values  of  the  piezoelectric 
coefficients.  The  higher  the  coefficient,  the  greater  the  induced  strain. 

For  the  AFC  lamina,  the  piezoelectric  coefficients  and  efei  are  given  as  material 
parameters.  For  the  G-1195  PZT  lamina,  a  single  constant  is  inadequate  to  relate  strain 
to  the  electric  field.  From  a  series  of  elastically-constrained  piezoceramics  tests,  it  was 
determined  that  the  value  of  ^33  depends  on  the  level  of  induced  strain  [8].  This  is  the 
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strain  that  appears  in  the  PZT  lamina  when  it  is  embedded  in  a  structure,  and  it  is  due  to 
the  piezoelectric  actuation.  Therefore,  a  secant  piezoelectric  coefficient  will  be  defined  as 

dh=  ±  (5-34) 

&  3 

where  Ai  is  the  actuation  strain  (that  physically  causes  induced  strains,  and  can  be  due 
to  piezoelectricity,  electrostriction,  or  even  thermal  expansion)  in  the  transverse  direction. 
An  approximate  linear  relation  can  be  set  up  between  the  induced  strain  and  the  applied 
electric  field,  using  Fig  A.l  from  reference  [10].  Another  linear  approximation  for  the 
appropriate  value  of  the  piezoelectric  constant  can  be  found  from  Figure  7,  of  reference 
[8].  The  value  of  the  transverse  piezoelectric  constant  d$i  can  be  approximated  as  10%  of 
g?33.  The  maximum  value  for  Ai  is  selected  as  the  maximum  strain  so  that  depoling  of  the 
PZT  fibers  is  avoided  [3]. 

Depoling  of  the  fibers  can  occur  if  the  applied  voltage  V  becomes  greater  than  the 
coercive  field  Ec.  During  the  manufacture  of  piezoceramics,  a  coercive  electric  field  is 
applied  across  the  fibers  to  align  the  PZT  crystals  into  an  initial  polarization.  If  the 
applied  voltage  during  operations  is  greater  than  this  electric  field  in  the  opposite  direction, 
depoling  of  the  fibers  can  take  place,  and  repoling  in  the  opposite  direction  will  occur.  If  the 
applied  field  is  aligned  with  the  initial  poling  direction,  depoling  will  not  take  place  even 
if  the  applied  electric  field  exceeds  the  coercive  field.  Therefore,  the  maximum  voltage 
applicable  to  the  AFC  lamina  will  be  ±1000  V  [5].  For  the  G-1195  PZT  lamina,  the 
applicable  voltage  per  thickness  is  ±750-1000  Vmm  [8], [15].  Considering  the  thickness  of 
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the  G-1195  PZT  lamina  listed  in  Table  (5.1),  the  maximum  voltage  that  can  be  applied  is 
±200V. 

Now  that  all  the  pieces  of  information  about  the  piezoelectric  coefficients  in  the 
material  principal  axes  are  given  for  both  PZT  laminae  through  Table  (5.1)  and  Eqn.  (5.34), 
using  the  appropriate  coordinate  transformations  we  find  the  PZT  strain  tensor  est  in  the 
structural  directions. 

5.3.2  PZT  Strain  Tensor  Transformation.  Figure  (5.5)  illustrates 
the  coordinate  rotation  from  the  material  principal  directions  into  the  structural  directions. 
The  principal  fiber  direction  of  the  PZT-composite  is  represented  by  the  principle  direction 
3,  and  the  structural  spanwise  direction  is  given  by  z.  The  principal  axes  3-1-2  are  rotated 


Figure  5.5  PZT  Strain  Tensor  Coordinate  Transformation 
by  angle  —9  about  the  2  axis  into  the  structural  axes.  Then,  the  principle-to-structural 
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transformation  (a  ’-2’  rotation)  matrix  C  is  given  by 

cos(9)  0  sin{9) 

C  =  0  10  (5-35) 

— sin{9)  0  cos(9) 

The  principle-to-structural  strain  tensor  transformation,  using  Eqns.  (5.33)  and  (5.35)  is 
given  by 

d3\  0  0 

est  ~  C  0  d3 i  0  C  E 3  (5-36) 

0  0  dss 

After  the  transformation,  we  obtain  the  PZT  strain  tensor,  as  well  as  the  PZT  strain  and 
shear  ep  and  in  the  structural  axes 

€x  0  £xz 

€$t  0  €y  0 

€zx  0  6Z 

where 

ez  -  ep  =  [d3\sin2{9)  +  d3zcos2(9 )]  E3 

exz  =  tzx  =  7 p  =  [d3\cos{9)sin{9)  -  d33sin(9)cos(9)]  E3 

are  to  be  used  in  Eqn.  (5.29). 


(5.37) 


(5.38) 
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5-4  Discussion 

5.4- 1  Composite-PZT  Beam  Torsion.  Evaluating  Eqns.  (5.9)  through 
(5.38)  —  simultaneously  for  the  given  geometry  and  composite-PZT  construction  —  is  a 
formidable  as  well  as  lengthy  task.  The  evaluation  of  the  full-,  and  the  partial-path  integrals 
could  easier  be  done  by  using  symbolic  solvers.  The  symbolic  solutions  of  Eqns.  (5.9) 
through  (5.29)  were  found  using  MathCad,  and  were  entered  into  a  Matlab  subroutine 
called  ANISOTORSIONJPZT  that  calculates  the  tip  twist  angle  of  the  non-homogeneous 
(see  Section  4.2)  anisotropic,  single-cell,  closed-section  box  beam  with  a  single  layer  of  PZT 
lamina  in  the  geometric  center  of  the  composite  laminate,  subjected  to  a  uniform  cross 
sectional-moment  about  the  structural  £  axis.  This  code  can  also  be  used  to  calculate  the 
tip  twist  of  a  homogeneous  anisotropic  beam  subjected  to  the  same  conditions,  since  in 
this  case  the  coupling  constant  @2  is  zero  (Section  4. 1.3. 2). 

The  Matlab  subroutine  called  LAYUP  of  Section  4.3  is  used  by  the  program  ANISO- 
TORSION_PZT  without  modifications. 

A  driver  program  called  COMPL_ANISO_PZT  was  written  in  order  to  define  the 
material  properties  in  the  1-2  principal  axes  of  the  composite  and  PZT  lamina  materials 
of  choice.  It  also  defines  the  single-cell  torquebox  geometry  obtained  from  the  program 
AREA  of  Chapter  II,  and  establishes  the  desired  loading  conditions.  Both  codes  are  listed 
in  Appendix  D. 

5.4- 2  PZT  Strain  Actuation.  The  Matlab  subroutine  called  ANISO- 
TORSIONJPZT_E  calculates  the  tip  twist  angle  of  the  non-homogeneous  (see  Section  4.2), 
anisotropic,  single-cell,  closed-section  box  beam  with  a  single  layer  of  PZT  lamina  in  the 
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geometric  center  of  the  composite  laminate,  subjected  to  a  uniform  electric  field  E:i  ap¬ 
plied  in  the  principal  PZT  fiber  direction  due  to  an  applied  voltage  V.  The  twist  angles  are 
calculated  for  a  range  of  fiber  orientation  angles,  input  by  the  user.  This  code  can  also  be 
used  to  calculate  the  tip  twist  of  a  homogeneous  anisotropic  beam  subjected  to  the  same 
conditions,  since  in  this  case  the  elastic  coupling  constant  is  zero  (see  Section  4. 1.3.2). 

The  program  ANISOTORSIONJPZT_E  uses  the  LAYUP  subroutine  of  Section  4.3 
without  modifications. 

The  transformation  of  the  PZT  strains  and  shears  in  the  principal  axes  into  the  strains 
and  shears  in  the  structural  axes  is  accomplished  by  the  subroutine  TRANSFORM.  Its 
inputs  are  the  PZT  principal  strain  tensor,  and  the  PZT  rotation  angle  with  respect  to  the 
structural  axes.  It  outputs  the  strain  tensor  in  the  structural  axes. 

A  driver  program  called  COMPL_ANISOJPZT_E  was  written  in  order  to  define  the 
material  properties  in  the  principal  axes  of  the  composite  and  PZT  lamina  materials  of 
choice.  It  also  defines  the  piezoelectric  properties  of  the  PZT  lamina  in  the  3-1  principal 
axes  and  establishes  the  desired  loading  conditions.  All  three  programs  and  subroutines 
are  listed  in  Appendix  E. 
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VI.  Results  and  Discussion 


6.1  Single and  Three- Cell,  Isotropic  Beam 

As  was  discussed  in  Chapter  III,  the  Matlab  code  SHEARFLOW3CELL  was  run 
in  order  to  gather  data  on  the  concentrated  torsional  moment  required,  and  the  resulting 
shear  flows  in  the  spars  and  skins  of  the  single-,  and  three-cell  closed  section,  isotropic 
(aluminum)  torque  box.  The  material  properties  of  the  2024-T3  aluminum  are  listed  in 
Table  (6.1).  The  code  was  run  with  the  original  geometric  properties  of  the  2-D  profile, 
obtained  from  the  Matlab  code  AREA  (see  Appendix  A)  for  angles  2  to  10  degrees  with 
2-degree  increments.  Interested  users  can  of  course  change  the  imposed  twist  angles  and 
investigate  the  relationship  between  the  angle  of  twist,  moment  and  shear  flows  required. 
For  the  given  angles,  Table  (6.2)  summarizes  the  results.  The  tabulated  values  of  the 
required  moments  versus  twist  angles  are  also  shown  (and  the  linear  relationship  can  be 
better  observed)  in  Figure  (6.1).  Comparing  the  moments  required  to  generate  the  angle 
of  twist  for  the  single-cell  beam  to  that  of  the  three-cell  beam  it  is  evident,  that  in  order  to 
generate  a  given  angular  twist,  a  larger  moment  is  required  for  the  three-cell  beam  than  for 
the  single-cell  beam.  This  is  due  to  the  increased  torsional  stiffness  (see  Eqn.  (3.29))  of  the 
multi-cell  beam,  compared  to  that  of  the  single-cell  beam.  This  is  the  reason  why  multi-cell 


Table  6.1  Isotropic  Material  Properties 


Aluminum  2024  T3 

Young’s  Modulus 
E  (Msi) 

Shear  Modulus 
G  (Msi) 

Poisson’s  Ratio 

V 

10.4 

3.86 

0.33 
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Table  6.2  Single-,  and  Three-Cell  Torque  Box  Results 


Single  Cell 

Three  Cell 

Twist  Angle 

Shearflow 

Moment 

ql 

q2 

q3 

Moment 

(deg) 

(lb/in) 

(ft-lb) 

(lb/in) 

(lb/in) 

(lb/in) 

(ft-lb) 

2.0 

170.6 

6517 

129.7 

165.9 

63.9 

8733 

4.0 

341.2 

13033 

259.3 

331.7 

17467 

6.0 

511.8 

19550 

388.9 

497.6 

26200 

8.0 

682.4 

26067 

518.6 

663.5 

255.5 

34934 

10.0 

853.0 

32583 

648.3 

829.4 

319.3 

43667 

sections  are  more  resistive  to  torsion.  Another  important  observation  from  Table  (6.2)  is 


Figure  6.1  Required  Moments  for  Single-,  and  Three-Cell  Beam  Torsion 


presented  in  Figure  (6.2).  Here,  the  shear  flows  generated  in  the  single-cell  section,  as  well 
as  the  shear  flows  generated  in  all  three  sections  of  the  three-cell  beam  are  plotted  against 
the  imposed  twist  angles.  It  can  be  seen  that  for  a  given  twist  angle  the  highest  skin  shear 
flow  is  generated  in  the  single-cell  beam.  This  is  due  to  the  fact  that  all  the  shear  is  carried 
by  this  one  section.  All  the  shear  flows  in  the  three-cell  beam  are  lower,  because  a  larger 
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Table  6.3  Simplified  Single- Cell  Torquebox  Results 


Single  Cell 

Twist  Angle 

Shearflow 

Moment 

(deg) 

(lb/in) 

(ft-lb) 

2.0 

140.3 

4404 

4.0 

280.5 

8808 

6.0 

420.7 

13213 

8.0 

561.0 

17617 

10.0 

701.2 

22021 

combined  area  of  three  cells  is  available  to  resist  torsion.  The  nose  and  the  tail  sections 
of  the  three-cell  beam  enclose  the  smallest  areas  (see  Table  (2.10)),  and  the  shear  flows 
generated  will  be  the  least  here.  The  shear  flow  in  the  midsection  of  the  three-cell  beam 
is  less,  but  almost  equal  to  that  of  the  single-cell  beam.  Because  the  enclosed  area  of  the 
mid  section  is  the  largest  of  the  three  cells,  it  will  take  up  most  of  the  shear  load  generated 
by  the  torsion. 


Figure  6.2  Shear  Flows  for  Single-,  and  Three-Cell  Beam  Torsion 
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To  check  the  accuracy  of  the  solution  derived  in  Section  4.1,  it  is  necessary  to  use  the 
simplified  torquebox  geometry  calculated  in  Section  2.5.2  for  obtaining  the  moments  and 
shearflows  for  a  given  twist  angle.  See  Table  (2.5)  for  the  dimensions  and  wall  thickness 
of  the  simplified  torquebox.  Running  the  Matlab  code  SHEARFLOW3CELL  (see  Ap¬ 
pendix  B)  and  concentrating  on  the  results  for  the  (now  simplified  dimension)  single-cell 
torquebox,  we  obtain  the  moments  and  shearflows,  listed  in  Table  (6.3). 

6.2  Generalized  Torsion  Solution 

6.2.1  Single- Cell,  Isotropic  Beam.  In  order  to  check  the  accuracy  of 
the  non-homogeneous,  anisotropic  solution  developed  in  Section  4.2  for  the  torsional  de¬ 
flection  of  the  single-cell  torquebox,  we  can  apply  it  to  the  isotropic,  homogeneous  case, 
developed  in  Section  4.1.  The  scope  of  the  analysis  covered  only  the  applied  concentrated 
cross-sectional  moment,  and  it  did  not  consider  the  other  distributed  forces  and  moments 
that  may  act  upon  the  section  (see  Figure  (4.1)).  As  was  discussed  in  Section  4. 1.3. 2,  in 
case  of  the  homogeneous  isotropic  beam,  the  coupling  elastic  constants  p2  =  0,  and  0:2  —  0. 
The  Matlab  program  COMPL_ANISO_COPMZ,  along  with  its  subroutine  ANISOTOR- 
SION.COMPZ  (see  Appendix  C)  was  run  for  the  isotropic  (aluminum)  case  by  setting  the 
material  properties  to  those  of  isotropic  aluminum  and  by  using  the  moments  required  to 
achieve  the  angle  of  twists  from  Table  (6.3).  The  results  were  then  tabulated  in  Table  (6.4) 
along  with  the  results  obtained  from  the  shearflow  solution  (see  Table  (6.3)). 

The  amount  of  twist  calculated  via  Libove’s  method  for  the  isotropic  single-cell  beam 
corresponded  well  to  the  angles  obtained  from  the  shearflow  solution  using  the  Matlab 
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Tab: 


ie  6.4  Libove’s  Method  for  Single  Cell  Isotropic  Beam 


Shearflow  Results 

Libove’s  Method 

Applied  Moment 
(ft-lb) 

Twist  Angle 
(deg) 

Twist  Angle 
(deg) 

2.0 

1.938  ! 

4.0 

3.875 

13213 

6.0 

5.813  j 

17617 

8.0 

10.0 

9.688 

code  SHEARFL0W3CELL  (Appendix  B).  While  the  trend  conserves  the  linearly  elastic 
assumption  made  in  Section  4.1.1  and  demonstrates  a  linearly  increasing  deflection  angle 
with  increasing  applied  moment,  it  is  only  3.1%  less  than  the  deflection  angles  set  for  the 
shearflow  solution. 


6.2.2  Single-Cell,  Anisotropic  Beam.  To  verify  the  accuracy  of  the 
fully  anisotropic  solution  of  Section  4.2  the  results  were  compared  to  the  linear  case  derived 
and  experimentally  tested  by  Romeo,  et.  al.  in  reference  [24].  The  authors  also  expanded 
the  theory  to  include  non-linear  twist  effects  due  to  the  non-linear  effective  shear  modulus 
of  the  skin  panels  [24].  They  used  a  graphite/epoxy,  single-cell,  rectangular  cross-section, 
composite  torquebox  under  pure  torsion.  The  M40/914,  [452/  —  452/O2/9O2].;  laminate  was 
applied  once  on  the  top  and  bottom  skins,  and  twice  on  the  main  and  rear  spars.  The 
geometric  and  material  properties  of  the  torquebox  are  listed  in  Table  (6.5).  This  geometry 
was  substituted  into  the  Matlab  code  COMPL_ANISO_COMPZ,  which  was  run  with  the 
subroutine  ANISOTORSION.COMPZ  for  several  values  of  the  applied  moments  listed  in 
reference  [24].  The  results  —  tabulated  in  Table  (6.6)  —  were  in  good  agreement  with  the 
results  published  by  Romeo. 
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Table  6.5  Composite  Material  Properties  and  Geometry 


Carbon/Epoxy  M40/914  Vicotex 

Ei  (Msi) 

E2  (Msi)  G\2  (Msi)  V12  L  (in)  H  (in)  tskin  (in) 

30.4 

1.0  0.6  0.31 

28.35  5.24  0.08 

Table  6.6  Composite  Beam  Torsion  Verification. 


Applied  Moment 
(ft-lb) 

Skin  Thickness 
(in) 

Romeo’s  Method 
Twist  Angle 
(deg) 

Libove’s  Solution 
Twist  Angle 
(deg) 

2950 

0.08 

0.210 

0.235 

3688 

0.275 

4425 

0.320 

0.353 

Having  verified  the  accuracy  of  the  torsion  solution  developed  in  Chapter  IV  by  ap¬ 
plying  it  to  both  the  isotropic  and  anisotropic  case,  the  program  COMPL_ANISO_COPMZ 
was  run  for  a  fully  anisotropic  case  using  the  original  baseline,  simplified  torquebox  dimen¬ 
sions  of  Section  2.5.2,  by  selecting  the  material  properties  of  the  graphite-epoxy  composite 
listed  in  Table  (6.7).  The  code  was  then  run  with  its  subroutines  LAYUP,  and  ANISOTOR- 
SION_COMPZ  in  order  to  calculate  the  engineering  properties  of  the  composite  laminate 
in  the  structural  axes,  to  determine  the  composite  compliance  matrix  S  in  the  structural 
axes,  and  to  calculate  the  angle  of  twist  of  the  anisotropic,  composite,  single-cell  beam. 
A  composite  laminate  of  Carbon/Epoxy  (AS4/3501-6)  with  arbitrary  lamina  fiber  orien¬ 
tation  angles  of  [0/  6  /45  /-Q  / 0]  was  chosen,  with  a  standard  lamina  thickness  of  0.005 
inch.  The  relevant  material  properties  of  the  lamina  are  included  in  Table  (6.7);  how¬ 
ever,  for  more  information  on  the  material,  please  refer  to  reference  [11].  By  varying  the 


Table  6.7  Composite  Material  Properties 


Carbon/Epoxy  AS4/3501-6 

Young’s  Modulus 

Shear  Modulus 

Poisson’s  Ratios 

Ei  (Msi)  |  E2  (Msi) 

G12  (Msi) 

vn 

f2i 

20.6  1.50 

1.04 

0.27 

0.02 
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Table  6.8  Composite  Beam  Torsion.  M— 4404  ft-lb. 


Composite 

Laminate 

Total  Thickness 
(in) 

Twist  Angle 
(deg) 

[0  15  45  -  15  0] 

0.025 

4.940 

[0  30  45  -  30  0] 

0.025 

3.651 

[0  45  45  -45  0] 

0.025 

3.513 

[0  60  45  -  60  0] 

0.025 

4.132 

[0  75  45  -  75  0] 

0.025 

4.950 

[0  90  45  -  90  0] 

0.025 

5.340 

composite  lamina  angle  6  from  0  to  90  degrees  with  2-degree  increments,  several  different 
fiber-orientation  laminates  were  run  to  investigate  how  the  different  laminae  angles  affect 
the  torquebox  twist  angles.  Table  (6.8)  summarizes  the  tip  twist  angles  of  several  selected 
laminate  construction  beams  subjected  to  M=4404  ft-lb  of  cross-sectional  moment.  Any 
other  distributed  or  concentrated  forces  and  moments  acting  on  the  beam  segment  were 
assumed  to  be  zero.  The  results  are  also  plotted  in  Figure  (6.3)  showing  the  familiar  pat¬ 
tern  of  the  change  in  the  angle  of  tip-twist  due  to  the  change  in  the  single  lamina  angle  9. 


Comparing  the  results  shown  in  Tables  (6.4)  and  (6.8)  we  come  to  the  conclusion 
that  the  model  composite  laminate  construction  torquebox  —  subjected  to  the  same  cross- 
sectional  moment  —  twists  more  than  the  isotropic  (aluminum)  beam.  This  is  due  to 
the  fact  that  while  the  graphite-epoxy  composite  laminate  has  greater  stiffness  than  the 
aluminum,  we  only  used  a  five-ply,  0.025  inch  thick  composite  laminate,  compared  to  the 
aluminum  beam  wall  thickness  of  0.072  inch,  0.049  inch,  and  0.048  inch  for  the  main 
spar,  rear  spar,  and  skin  thickness  respectively,  which  resulted  in  lower  torsional  stiffness 
compared  to  the  aluminum  beam.  The  composite  wall  thickness  can  be  increased  by 
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increasing  the  number  of  laminae  or  the  laminates  up  to  the  point  where  it  satisfies  the 
design  stress  requirements  without  a  single  lamina  failure. 


Figure  6.3  Twist  Angles  for  the  All-Composite  Torquebox 

6.2.3  Single-Cell ,  Anisotropic  Composite-PZT  Beam.  The  Mat- 
lab  code  COMPL_ANISO_PZT,  developed  in  Section  5.1  to  calculate  the  twist  angle  of 
the  anisotropic  composite  beam,  was  modified  so  that  a  layer  of  PZT  embedded  in  the 
geometric  center  (mid-plane)  of  the  laminate  can  be  accounted  for  while  calculating  the 
engineering  properties  of  the  laminate  in  the  structural  axes.  This  was  achieved  by  or¬ 
dering  the  laminate  material  properties  in  the  material  principal  directions  into  a  vector, 
containing  the  PZT  material  properties  (E,G }u)  as  an  element.  Because  the  PZT  layer 
was  assumed  to  be  embedded  in  the  symmetric  center  of  the  laminate,  the  PZT  material 
properties  were  placed  in  the  middle  of  the  vector.  The  two  PZT  actuator  laminae  consid- 
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ered  were  introduced  in  Section  5.1.1.  The  mechanical  properties  of  the  AFC  lamina  —  in 
form  of  compliance  and  plane-stress  stiffness  matrix  —  were  obtained  from  a  combination 
of  experimental  stress-strain,  and  clamped-actuation  testing  by  Bent,  et.  al.  [5].  These 
plane-stress  compliance  values  were  used  to  calculate  the  material  properties  of  the  AFC 
lamina  in  the  principal  directions.  The  material  properties  of  the  G-1195  PZT  lamina  were 
obtained  from  reference  [3],  [8],  and  [10].  The  engineering  properties  of  both  PZT  laminae 
in  the  structural  directions  were  again  obtained  via  the  LAYUP  subroutine.  For  more 
details,  please  refer  to  the  Matlab  code  COMPLJ\.NISO_PZT  included  in  Appendix  D. 

According  to  the  argument  made  in  Section  5.1.1,  the  PZT  lamina  was  assumed  to 
be  at  a  45  degree  angle  with  respect  to  the  structural  z  axis  (see  Figure  (5.2)),  simply 
replacing  one  layer  of  the  graphite/epoxy  lamina.  Because  the  engineering  properties  of 
the  composite-PZT  laminate  are  different  —  in  the  material  1-2  axes,  as  well  as  in  the 
structural  z-x  axes  —  than  those  of  the  graphite-epoxy  composite  laminate  alone,  it  was 
necessary  to  calculate  all  the  elastic  coupling  coefficients  separately  for  the  composite  and 
composite-PZT  hybrid  laminate.  Therefore,  in  addition  to  the  compliance  coefficients 
a\.  «2,  and  0:3  of  the  composite  laminate,  we  will  have  ap  1,  ap 2,  and  apz  for  the  hybrid 
composite-PZT  laminate.  Similarly,  we  obtain  the  values  of  the  elastic  coupling  coefficients 
jS pi,  /3P2,  and  /3P3  for  the  composite-PZT  laminate.  These  coefficients  had  to  be  substituted 
in  the  appropriate  locations  in  all  the  partial  and  full-path  integrals  detailed  in  Chapter  IV. 
These  integrals  were  solved  symbolically  using  MathCad,  and  the  results  were  transferred 
to  the  Matlab  subroutine  ANISOTORSION_PZT.  This  subroutine,  when  called  by  the 
driver  program  COMPL_ANISO_PZT,  calculates  the  tip  twist  angle  due  to  a  cross-sectional 
moment  of  M=4404  ft-lb  of  the  single-cell,  closed-section,  anisotropic  torquebox  with  PZT 
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Table  6.9  PZT-Composite  Beam  Torsion.  M— 4404  ft-lb. 


AFC  Composite 

G-1195  PZT 

Composite 

Total  Thickness 

Twist  Angle 

Total  Thickness 

Twist  Angle 

Laminate 

(in) 

(deg) 

(in) 

(deg) 

[0  PZT  0] 

0.0165 

18.3149 

0.015 

19.8341 

[0  15  PZT  -  15  0] 

0.0265 

6.3218 

0.025 

6.6210 

[0  30  PZT  -  30  0] 

0.0265 

0.025 

6.8166 

[0  45  PZT  -  45  0] 

0.0265 

6.6154 

0.025 

6.6587 

[0  60  PZT  -  60  0] 

0.0265 

8.2119 

0.025 

8.3173 

[0  75  PZT  -  75  0] 

0.0265 

10.4376 

0.025 

10.7839 

[0  90  PZT  -  90  0] 

0.0265 

11.5805 

0.025 

12.1206 

actuator  lamina  embedded  in  the  center  layer  of  the  top  and  bottom  skin  graphite-epoxy 
composite  laminate.  The  program  was  run  for  the  same  composite  laminate  construction 
detailed  in  Table  (6.8)  with  the  center  45-degree  graphite-epoxy  layer  substituted  with 
one  layer  of  PZT  lamina  of  the  same  orientation.  Again,  the  two  PZT-composite  laminae 
used  for  comparisons  were  detailed  in  Section  5.1.1.  The  fiber  orientation  angle  6  of  the 
graphite-epoxy  lamina  was  varied  as  before  (see  Table  (6.8)),  and  the  results  are  tabulated 
in  Table  (6.9)  and  plotted  in  Figure  (6.4). 

Both  the  AFC  and  the  G-1195  PZT-Composite  beam  demonstrated  comparable  twist 
angles  for  identical  cross-sectional  moment  for  all  values  of  the  changing  lamina  angle. 
However,  both  beams  —  with  the  PZT  layer  embedded  —  exhibited  greater  tip  twist  angles 
for  a  given  substrate  lamina  orientation  angle  than  the  composite-alone  construction  beam. 
This  leads  to  the  conclusion  that  both  PZT  actuator  laminae,  when  embedded  in  the  center 
of  the  composite  laminate,  lowered  the  beam’s  torsional  stiffness,  compared  to  that  of  the 
composite-alone  laminate,  resulting  in  greater  beam  tip  twist  angles.  The  reduction  of 
beam  stiffness  was  expected  since  the  stiffness  of  both  actuator  laminae  are  one  fourth  of 
the  stiffness  of  the  substrate  graphite/epoxy,  and  the  shear  modulus  of  both  laminae  are 
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Figure  6.4  Twist  Angles  for  the  PZT-Composite  Torqueboxes 


an  order  of  magnitude  less  than  that  of  the  substrate  composite.  While  only  static  torsion 
is  considered,  the  reduction  in  torsional  stiffness  may  very  well  be  a  desirable  effect,  since 
a  lower  electric  field  is  necessary  to  generate  the  desired  angles.  When  dynamic  effects  are 
also  considered,  lowering  the  beam’s  torsional  stiffness  could  significantly  lower  the  wing’s 
flutter  velocity. 

6.2.4  Single- Cell,  Anisotropic ,  Composite-PZT  Beam  with  Strain 
Actuation.  The  driver  code  COMPL_ANISO_PZT  was  modified  in  order  to  account 
for  the  piezoelectric  actuation  of  the  embedded  PZT  lamina.  The  actuation  was  achieved 
by  applying  an  electric  field  in  the  PZT  lamina  poling  direction  (see  Figure  (5.3),  or  Fig¬ 
ure  (5.4)).  The  piezoelectric  strain  and  shear  generated  by  the  electric  field  was  accounted 
for  by  extending  the  constitutive  relations  and  re-deriving  Libove’s  method  in  Section  5.2. 
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Tal 


ble  6.10  Strain  Actuation  Twist  Angles  -  AFC  Laminate 


Composite 

Laminate 

100 

Appli 

250 

ed  Voltag 
500 

;e  (V) 

750 

1000 

[0  15  45  -  15  0] 

0.0299 

0.0748 

0.1495 

0.2243 

0.2991 

[0  30  45  -  30  0] 

0.0270 

0.0674 

0.1348 

0.2022 

0.2696 

[0  45  45  -  45  0] 

0.0282 

0.0704 

0.1409 

0.2113 

0.2818 

[0  60  45  -60  0] 

0.0323 

0.0807 

0.1615 

0.2422 

0.3229 

[0  75  45  -  75  0] 

0.0375 

0.0936 

0.1873 

0.2809 

0.3746 

[0  90  45  -  90  0] 

0.0400 

0.0999 

0.1998 

0.2997 

0.3996 

The  solution  was  coded  as  the  Matlab  subroutine  ANISOTORSION_PZT_E.  The  piezo¬ 
electric  strain  tensor  in  the  material  principal  3-1  axes  was  transformed  to  the  structural 
z-x  axes  by  the  Matlab  subroutine  TRANSFORM,  using  the  methods  of  Section  5.3. 

The  Matlab  code  COMPL_ANISO_PZT_E  was  run  for  the  AFC  and  G-1195  PZT 
laminate  cases,  using  several  values  of  the  applied  electric  field  E%,  and  substrate  lamina 
orientation  angles.  The  results  for  the  strain  actuation  of  the  AFC  laminate  beam  are 
tabulated  in  Table  (6.10),  and  —  using  the  variable  lamina  angle  6  as  the  parameter 
—  the  tip  twist  angles  were  plotted  against  the  applied  voltage  in  Figure  (6.5).  The 
linear  trend  is  conserved,  and  the  effect  of  the  lamina  angle  on  the  tip  twist  for  a  given 
applied  voltage  is  conveniently  observed.  For  any  given  lamina  angle  9,  an  increase  in  the 
applied  voltage  increases  the  beam  tip  twist  angle.  For  a  given  applied  voltage,  the  tip 
twist  angle  decreases  until  the  lamina  angle  reaches  30  degrees.  Further  increase  in  the 
lamina  angle  increases  the  tip  twist,  indicating  that  the  host  structure’s  torsional  stiffness 
is  reduced  when  the  lamina  angle  is  greater  than  30  degrees.  To  better  observe  the  effect  of 
the  variable  angle  lamina,  and  to  more  accurately  determine  the  variable  substrate  angle 
where  the  minimum  twist  occurs,  the  tip  twist  angles  of  the  AFC  laminate  torquebox  were 
also  plotted  against  the  variable  lamina  angle  9  in  Figure  (6.6),  using  the  applied  voltage 
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Figure  6.5  Twist  Angles  Due  to  Strain  Actuation  -  AFC  Laminate 


as  the  parameter.  The  plot  demonstrates  how  the  variable  lamina  angle  affects  the  tip 
twist  for  a  given  value  of  applied  voltage.  By  inspection,  the  minimum  twist  for  any  given 
applied  voltage  occurs  at  approximately  0—33  degrees.  Any  increase  or  decrease  in  the 
lamina  orientation  angle  will  result  in  an  increase  in  the  tip  twist  angle. 

The  results  for  tip  torsion,  using  the  G-1195  PZT-composite  hybrid  laminate,  are 
tabulated  in  Table  (6.11).  The  tip  twist  angles  are  greater  than  those  of  the  AFC  laminate 
beam,  and  they  correspond  well  to  the  trend  observed  in  Table  (6.10).  This  is  partly  at¬ 
tributed  to  the  piezoelectric  constants  of  the  G-1195  lamina  which  are  significantly  greater 
than  those  of  the  AFC  lamina  (see  Table  (5.1)).  The  tip  twist  angles  of  the  G-1195  PZT 
laminate  torquebox  are  again  plotted  against  the  applied  voltage  in  Figure  (6.7),  using 
the  variable  substrate  lamina  angle  0  as  the  parameter.  The  linear  trend  is  not  conserved 
as  before,  because  the  induced  strain  increase  with  applied  electric  field,  which  in  turn 
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Figure  6.6  Twist  Angles  Due  to  Strain  Actuation  -  AFC  Laminate 


Table  6.11  Strain  Actuation  Twist  Angles  -  G-1195  PZT  Laminate 


Composite 

Laminate 

50 

Appli 

75 

ed  Voltaj 
100 

5e  (V) 

150 

200 

[0  15  45  -  15  0] 

0.1184 

0.1899 

0.2695 

0.4533 

0.6698 

[0  30  45  -  30  0] 

0.0992 

0.1590 

0.2257 

0.3796 

0.5609 

[0  45  45  -  45  0] 

0.1071 

0.1718 

0.2438 

0.4101 

0.6060 

[0  60  45  -  60  0] 

0.1340 

0.2149 

0.3050 

0.5130 

0.7579 

[0  75  45  -  75  0] 

0.1676 

0.2688 

0.3816 

0.6418 

0.9484 

[0  90  45  -90  0] 

0.1840 

0.2951 

0.4188 

0.7044 

1.0409 

increase  the  value  of  the  piezoelectric  constant  in  the  poling  (^33)  direction.  For  any  given 
substrate  lamina  angle,  the  tip  twist  angle  increases  with  increasing  applied  voltage,  while 


for  a  given  applied  voltage,  the  tip  twist  angle  decreases  until  the  lamina  angle  reaches  30 
degrees.  Similar  to  the  AFC  laminate  case,  further  increase  in  the  lamina  angle  increases 
the  tip  twist,  indicating  that  the  host  structure’s  torsional  stiffness  is  reduced  when  the 
composite  lamina  angle  is  greater  than  30  degrees. 
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Figure  6.7  Twist  Angles  Due  to  Strain  Actuation  -  G-1195  PZT 


The  tip  twist  angles  of  the  G-1195  PZT  laminate  torquebox  were  also  plotted  against 
the  variable  lamina  angle  6  in  Figure  (6.8),  using  the  applied  voltage  as  the  parameter. 
The  plot  demonstrates  how  the  variable  lamina  angle  affects  the  tip  twist  for  a  given  value 
of  applied  voltage.  By  inspection,  the  minimum  twist  for  any  given  applied  voltage  occurs 
at  approximately  9= 33  degrees.  Any  increase  or  decrease  in  the  lamina  orientation  angle 
will  result  in  the  increase  of  the  tip  twist  angle.  The  recorded  tip  twist  angles  for  a  given 
applied  voltage  again  proved  to  be  greater  for  the  composite  lamina  orientation  angle  of 
90  degrees,  than  at  0  degrees.  The  same  result  was  observed  when  using  the  AFC  actuator 
lamina. 
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Figure  6.8  Twist  Angles  Due  to  Strain  Actuation  -  G-1195  PZT 


6.3  Discussion 

The  results  presented  for  the  single-,  as  well  as  the  three-cell  isotropic  beam  well 
demonstrated  the  effects  of  the  higher  torsional  stiffness  afforded  by  the  multiple  cell 
section,  resulting  in  less  tip  twist.  When  less  torsional  displacement  is  desired,  the  multiple¬ 
cell  construction  offers  better  results;  however,  the  enclosed  area  will  be  greater,  requiring 
a  larger  structure.  The  solutions  to  the  isotropic,  single-cell  trapezoid  torquebox  (see 
Table  (6.3))  were  prepared  to  verify  the  accuracy  of  the  anisotropic  torsion  model,  when 
applied  to  full  isotropy. 

The  solution  to  the  torsion  of  the  single-cell,  trapezoid  cross-section,  anisotropic 
beam  —  simplified  to  isotropic  material  properties  —  matched  the  results  obtained  from 
the  isotropic  solution  within  3.1%.  The  results  proved  the  validity  of  the  torsion  model, 
when  used  with  a  uniform  cross-sectional  moment  acting  on  the  section.  The  full  aniso- 
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tropic  model  —  with  cross  section  and  other  dimensions  modified  to  match  the  torquebox 
dimensions  of  reference  [24]  —  predicted  torsional  displacements  close  to  those  calculated 
by  the  linear  theory  presented  by  reference  [24]  (see  Table  (6.6)).  The  anisotropic  torsion 
model  predicted  higher  angular  tip  displacements  than  those  listed  in  reference  [24].  The 
difference  can  be  attributed  to  the  approximation  used  in  reference  [24]  by  extending  the 
Bredt-Batho  theorem  of  isotropic  beams  to  approximate  the  twist  of  anisotropic  beams. 
Both  linear  models,  however,  underpredict  the  experimentally  determined  tip  twist  angles 
of  the  anisotropic  beam  due  to  assuming  linear  strains  and  rotations  (Eqns.  (4.20)  and 
(5.5)),  as  well  as  constant  shear  stiffness.  The  accuracy  of  the  method  of  Romeo,  et.  al. 
will  increase  when  the  incomplete  diagonal  shear  stress  field  in  panels  operating  in  the 
post-buckling  phase  is  accounted  for  [24].  While  the  model  of  anisotropic,  single-cell  beam 
torsion  —  derived  in  Chapter  IV  —  offers  higher  fidelity  in  the  linear  regime,  its  accuracy 
can  further  be  increased  by  accounting  for  non-linear  translational  strains,  as  well  as  non¬ 
linear  curvatures  about  the  structural  axes.  The  accuracy  can  be  even  further  increased  by 
accounting  for  warping  functions  of  various  degrees  (linear,  or  non-linear),  as  was  discussed 
in  Section  4.1. 3.1. 

The  PZT  lamina  was  embedded  in  the  host  structure  (substrate),  and  the  engineering 
properties  of  the  hybrid  laminate  were  recalculated  so  that  the  modified  compliance  values 
could  be  used  to  calculate  the  material  elastic  coupling  coefficients.  Both  PZT  laminae 
(the  AFC  of  reference  [5],  and  the  G-1195  PZT  lamina  of  reference  [3])  were  accounted 
for  by  using  their  respective  material  elastic  and  piezoelectric  properties.  The  torquebox 
tip-torsion  angles  due  to  the  applied  voltage  were  recorded,  tabulated,  and  presented  in 
a  parametric  plot,  using  first  the  variable  composite  ply  angle,  then  the  applied  voltage 
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as  a  parameter.  The  strain  actuation  plots  of  the  AFC  PZT  lamina,  using  the  composite 
substrate  lamina  angle  as  parameter  demonstrated  a  linear  relation  between  applied  voltage 
and  twist  angle.  The  strain  actuation  plots  of  the  G-1195  PZT  lamina,  using  the  composite 
substrate  lamina  angle  as  a  parameter,  demonstrated  a  non-linear  relation  between  applied 
voltage  and  twist  angle.  This  was  due  to  the  piezoelectric  constant,  which  is  a  function 
of  the  induced  strain  of  the  PZT  lamina.  The  results  probably  still  underpredicted  any 
experimental  results  due  to  the  arguments  made  in  the  previous  paragraph.  Also,  the 
angles  obtained  due  to  strain  actuation  are  most  probably  insufficient  to  provide  adequate 
aerodynamic  control  of  the  UAV  designed  in  Chapter  II  at  any  airspeed  within  the  flight 
envelope.  The  roll  authority  provided  by  the  above  strain  actuation  results  can  be  evaluated 
by  any  available  theory,  such  as  lifting-line  theory,  or  numerical  (panel)  methods. 

The  variation  of  twist  angle  as  a  function  of  the  composite  substrate  lamina  angle  — 
using  the  applied  voltage  as  a  parameter  —  produced  patterns  similar  to  those  obtained 
from  the  torsion  of  the  pure  composite  beam,  and  the  torsion  of  the  PZT-composite  beam 
Sections  6.2.2  and  6.2.3.  The  minimum  twist  angle  due  to  constant  voltage  strain  actu¬ 
ation  of  both  PZT  laminated  torqueboxes  occurred  at  the  composite  substrate  angle  of 
approximately  33  degrees. 

Comparing  the  results  tabulated  in  Tables  (6.11)  and  (6.10),  and  plotted  in  Fig¬ 
ures  (6.6)  and  (6.8),  the  torquebox  using  the  G-1195  PZT  lamina  achieved  higher  twist 
angles  with  lower  applied  voltages  than  the  torquebox  equipped  with  the  AFC  lamina. 
Though  both  piezoelectric  constants  of  the  G-1195  lamina  are  positive  (creating  strains  of 
equal  signs  for  a  given  applied  electric  field),  they  are  considerably  greater  than  those  of 
the  AFC  lamina.  This  confirms  the  statement  made  in  Section  5.1,  that  the  PZT  lamina  of 
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high  e?33  will  require  less  voltage  to  produce  the  necessary  strains.  The  results  also  validate 
the  claim  that  the  higher  the  stiffness  the  larger  the  fraction  of  the  strain  produced  by  the 
electric  field  that  is  transferred  to  the  host  structure.  That  is,  the  actuator  lamina  will  be 
compressed  less,  while  the  substrate  will  be  strained  more. 
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VII.  Conclusions  and  Recommendations 


7. 1  Summary 

This  thesis  detailed  a  method  of  analyzing  and  modeling  a  fully  anisotropic,  single- 
cell,  closed-section,  generic  trapezoid  cross-section  torquebox  under  torsional  deformation, 
using  a  full  complement  of  air  loads.  The  conceptual  design  of  a  proposed  low-observable, 
medium-range  and  endurance  reconnaissance  UAV  was  accomplished  in  order  to  develop 
the  baseline  geometry  for  the  three-cell,  as  well  as  the  single-cell,  isotropic  wing  torquebox, 
as  well  as  to  provide  prototype  dimensions  for  the  RCS  comparative  study  in  reference 
[25].  The  simplified  single-cell  torquebox  dimensions  were  generated  by  assuming  equal 
length,  flat-surface  wing  skins  between  the  spars  as  to  facilitate  the  derivation  of  the  fully 
anisotropic  torsion  model.  The  single-cell  beam  dimensions  were  then  used  to  determine 
the  uniform  twisting  moment  required  to  generate  the  given  amount  of  tip  twist  angles  of 
the  isotropic  torquebox,  as  predicted  by  the  Bredt-Batho  theorem. 

A  torsion  model  was  then  derived  for  both  the  homogeneous  and  non-homogeneous 
anisotropic  structures.  The  fully  anisotropic  case  was  developed  based  on  modifying  Li- 
bove’s  method  using  a  thin- walled,  linearly  elastic,  fully  anisotropic,  trapezoid  cross-section 
beam.  The  accuracy  of  the  torsion  solution  using  a  fully  isotropic  case  was  compared  to 
the  solution  obtained  for  isotropic  beam  torsion  using  the  Bredt-Batho  theorem.  The 
fully  anisotropic  case  (modified  to  a  rectangular,  non-constant  thickness  cross-section  box 
beam)  was  compared  to  the  results  published  by  Romeo,  et.  al.  for  anisotropic  laminate 
box  beams. 
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The  investigation  extended  the  analysis  to  include  a  PZT-composite  actuator  lamina 
in  the  mid-plane  of  the  composite  host  structure’s  skin,  and  to  model  the  piezoceramic 
actuation  through  applied  piezoceramic  strain  and  shear.  The  effect  of  the  PZT-composite 
lamina  on  the  torsional  stiffness  of  the  torquebox  was  demonstrated  by  applying  a  uniform 
cross-sectional  moment  to  the  hybrid  composite  torquebox. 

The  modeling  of  the  piezoceramic  actuation  was  accomplished  by  extending  Libove’s 
method  derived  for  the  anisotropic  composite  torquebox,  to  include  piezoelectric  strain  and 
shear  in  the  strain-stress  constitutive  relations.  The  torsional  model  was  then  re-derived 
to  include  the  effect  of  PZT  strain  and  shear  present  in  the  top  and  bottom  skin  surfaces. 
The  PZT  embedded  surfaces  were  strain  actuated  by  applying  an  electric  field  to  the  PZT 
actuators,  inducing  torsion  of  the  host  structure.  The  voltage  generated  tip  torsion  of  the 
torquebox  was  verified  by  recording  the  angles  of  twist  due  to  a  range  of  applied  voltages. 
The  veracity  of  the  solution  was  demonstrated  by  varying  the  lamina  angle  of  two  layers 
of  the  graphite-epoxy  composite,  and  the  resulting  tip  torsion  angles  were  recorded. 

The  trend  of  twist  angles  versus  substrate  lamina  angles  corresponded  to  the  trend 
observed  for  the  composite  beam  torsion  due  to  applied  moment.  The  torsion  of  the 
torquebox  due  to  strain  actuation  using  two  different  PZT  laminae  were  compared. 

7.2  Conclusions 

This  research  has  covered  a  span  of  over  six  months,  starting  from  the  conceptual 
design  and  iteration  of  the  proposed  UAV.  The  shearflow  solution  to  the  closed-section, 
single-cell  and  three-cell  isotropic  torquebox  torsion  was  developed  by  using  the  easily 
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applicable  Bredt-Batho  theorem.  The  derivation  of  the  anisotropic  torsion  model,  however, 
was  achieved  only  after  overcoming  considerable  conceptual  as  well  as  algebraic  difficulties 
in  establishing  the  appropriate  mathematical  boundary  conditions  and  correct  limits  of 
integration  for  the  simplified  single-cell  torquebox.  Some  difficulty  was  encountered  in 
arriving  at  the  correct  plane-stress  constitutive  elements  that  were  used  to  model  the  strain 
and  shear  of  the  composite  laminate.  The  generating  of  the  correct  symbolic  solutions  of 
over  30  partial,  as  well  as  the  full-path  integrals  for  the  piecewise  continuous  surfaces 
—  that  were  indispensable  to  develop  the  code  that  automated  the  solution  —  was  a 
particularly  slow  and  tedious  process,  after  which  the  nature  of  the  research  process  became 
painfully  obvious: 

The  research  process  is  extremely  non-linear  if  not  discontinuous.  A  year’s  or 

decade’s  worth  of  work  can  pay  off  in  one  day.  ...  Be  patient  and  persistent  [6] 

Once  the  process  was  automated  via  the  Matlab  routines,  the  numerical  results  for  the 
different  cases  became  easier  to  generate. 

Nevertheless,  the  procedure  had  to  be  completely  repeated  for  the  PZT-composite 
laminate  case,  since  the  plane-stress  constitutive  equations  were  modified  to  account  for 
the  PZT  lamina,  as  well  as  for  the  PZT  strain  and  shear.  Once  the  PZT-composite  solution 
to  the  integrals  were  completed  using  the  appropriate  elastic  coupling  coefficients  of  the 
composite,  and  that  of  the  PZT-composite  hybrid  laminate,  the  existing  Matlab  routines 
had  to  be  modified  accordingly  in  order  to  account  for  the  respective  changes  in  the  path 
integrals.  Only  then  it  became  routine  to  generate,  run,  and  record  sample  cases  to  test 
the  applicability  of  the  torsion  model  of  the  anisotropic  PZT  actuated  composite  laminate. 
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The  model  was  developed  in  a  way  that  is  easily  modified  to  include  different  com¬ 
posite  materials,  any  different  PZT  actuator  laminas,  as  well  as  any  trapezoid  cross-section 
areas  (that  can  include  rectangles,  squares,  or  even  to  approximate  plates),  thereby  pro¬ 
viding  an  extremely  fast  and  convenient  tool  for  initial  theoretical  prediction  of  torquebox 
torsion  in  the  linear  regime. 

7. 3  Recommendations 

The  thesis  achieved  the  goals  set  out  in  Section  1.3;  however,  it  could  be  further 
expanded  by  investigating  several  other  issues  that  are  relevant  to  the  topic  of  beam  torsion. 
The  study  can  be  continued  by  establishing  the  shear-moment  diagram  of  any  desired 
flight  condition  for  any  of  the  UAV  designs  of  Chapter  II,  and  the  incremental  shear  forces 
and  bending  moments  along  the  span  can  be  imposed  upon  the  single-cell  isotropic,  or 
anisotropic  torquebox.  Strain  actuation  can  be  added  to  the  structure  if  so  desired.  The 
angle  of  twist  for  a  given  particular,  or  range  of  loading  condition  can  be  calculated. 

Also,  the  torqueboxes  can  be  designed  to  the  given  maximum  flight  loads,  so  that 
actual  composite  skin  and  spar  thickness,  along  with  their  respective  lamina  orientation 
can  be  used  for  calculating  the  constitutive  relations,  as  well  as  the  angles  of  twist.  The 
design  methods  of  Section  2.4  will  have  to  be  abandoned,  as  they  only  apply  to  isotropic 
materials. 

The  single-cell  torquebox  solution  can  also  be  expanded  by  considering  a  three-cell 
(or  other  multi-cell)  beam.  This  can  be  achieved  by  supplementing  the  equations  of  static 
equivalence  of  Section  4. 1.3. 2  by  the  compatibility  equations  that  require  equal  amount 


7-4 


of  rotations  for  all  sections.  The  PZT  elements  then  can  be  applied  to  the  surfaces  of 
interest,  and  a  full  complement  of  shear  loads,  bending  moments  and  strain  actuation  can 
be  assumed  along  the  span  by  using  the  methods  of  Chapters  IV  and  V. 

After  obtaining  the  solution  to  the  multi-cell  torquebox,  the  six  UAV  designs  devel¬ 
oped  can  be  used  to  conduct  trade  studies  for  PZT  actuation  requirements  for  different 
takeoff  weights,  aspect  ratios,  etc.  Relations  can  be  determined  between  maximum  gross 
takeoff  weight  and  strain  actuation  power  requirements,  given  certain  design  variables  (such 
as  AR,  wingspan,  etc.)  are  held  constant. 

Once  the  trade  studies  are  completed,  the  results  could  be  extrapolated  within  a 
reasonable  range  (to  be  determined)  and  the  entire  process  could  be  automated,  so  that 
intermediate  values  of  the  independent  variable  (weight,  AR,  wingspan,  etc.)  can  be  linked 
to  the  dependent  variables  (actuation  strain,  power  required,  etc.). 

Finally,  the  author  recommends  an  experimental  verification  of  the  results  obtained, 
by  constructing  two  trapezoid  torqueboxes  of  any  desired  dimension  (although  I  recommend 
the  two  to  be  the  same  in  dimensions  to  ease  the  burden  as  well  as  lessen  the  expense  of 
manufacturing),  one  specimen  with,  the  other  without  the  PZT-composite  actuator.  The 
experiment  would  then  draw  conclusions  about  the  accuracy  of  the  results  presented  in 
Chapter  VI,  as  well  as  make  suggestions  for  possible  corrections  and  improvements. 


7-5 


Appendix  A.  Airfoil  Profile  Parameters 


Codes 


function [Al, 11, A2, 121, 122, A3, 131, 132, hi, h2] =area(cll ,c22,c) ; 

'/.Written  by  Capt  Peter  Cseke,  Jr.  -  Fall  1999  -  AFIT/ENY 
7, Calculates  the  geometric  properties  of  the  airfoil  profile 
'/.Inputs:  NLF0215.m  subroutine  for  airfoil  input  data 

'/,  interp.m  subroutine  for  interpolating  values 

'/.Outputs:  Geometric  properties  of  airfoil  (areas,  skin  lengths,  spar  heights) 

'/.Call  the  profile  function  for  non-dimensional  data 

[xcupper , zcupper , xclower , zclower] =NLF0215 (xcupper , zcupper , xclower , zclower)  ; 

’/.Count  the  length  of  the  vectors  for  the  x  coordinates 
'/,(=length  of  y  coordinates) 
nupper=length (xcupper) ; 
nlower=length (xclower) ; 

croot=c; 
cl=cll*croot ; 
c2=c22*croot ; 

'/.Calculate  the  dimensional  parameters  from  chord  and  profile  data 
xcupper=xcupper . *croot ;  zcupper=zcupper . *croot ; 
xclower=xclower . *croot ;  zclower=zclower . *croot ; 

7#*  *************************************************************** 

'/.Calculate  the  nose  cone  upper  area  (dAlup)  ,  and  upper  skin  length  (dllup) 


dAlup=0 .0; 
dllup=0.0; 

while  xcupper(u)<=cl 


xl=xcupper (u) ; 
x2=xcupper (u+1) ; 
yl=abs (zcupper (u) ) ; 
y2=abs (zcupper (u+1) ) ; 
x=xl+(x2-xl) / 2 ; 

[y] =interp(xl ,x2 ,x,yl ,y2) ; 

if  x2>cl 

dAlup=dAlup+ (cl-xl) *y ; 
dllup=dllup+sqrt ( (cl-xl) "2+ (y-yl) "2) ; 
else 

dAlup=dAlup+ (x2-xl) *y ; 
dllup=dllup+sqrt ( (x2-xl) "2+ (y-yl) "2) ; 

end 

u=u+l ; 

end 

'/.Reset  counter  so  that  first  index  starts  where  we  left  off 
u=u-l ; 

dAlup;  7, Display  result  if  so  desired 
dllup;  '/.Display  result  if  so  desired 

'/.Calculate  the  nose  cone  lower  area  (dAllow)  and  lower  skin  length  (dlllow) 

j=1i 

dAllow=0 . 0 ; 
dlllow=0. 0; 
while  xclower ( j) <=cl 
xl=xclower (j) ; 
x2=xclower ( j+1) ; 
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yl=abs (zclower ( j ) ) ; 
y2=abs(zclower(j  +  l) ) ; 
x=xl+(x2-xl)/2; 

[y]=interp(xl,x2,x,yl ,y2) ; 

if  x2>cl 

dAllow=dAllow+(cl-xl)*abs(y) ; 
dlllow=dlllow+sqrt ( (cl-xl) ~2+(y-yl) “2) ; 
else 

dlllow=dlllow+sqrt ( (x2-xl) "2+ (y2-yl) "2) ; 
dAllow=dAllow+ (x2-xl) *abs (y) ; 

end 

j-j+i; 

end 

j=j-l;  '/, Reset  counter  so  that  first  index  starts  where  we  left  off 

xclower(j);  '/.Display  result  if  so  desired 

dAllow;  '/.Display  result  if  so  desired 

dlllow;  '/.Display  result  if  so  desired 

'/.Calculate  the  sum  of  the  nose  cone  areas  (dAl)  and  skin  lengths  (11) 

Al=dAlup+dAllow 

ll=dllup+dlllow 

’/, Calculate  main  torque  box  upper  area  (dA2up)  and  upper  skin  length  (121) 
dA2up=0 . 0 ; 

121=0.0; 

while  xcupper (u) <=c2 
xl=xcupper (u) ; 
x2=xcupper (u+1) ; 
yl=zcupper (u) ; 
y2=zcupper (u+1) ; 
x=xl+ (x2-xl) /2 ; 
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[y]=interp(xl ,x2,x,yl ,y2) ; 


if  x2>c2  '/.If  beyond  spar,  go  back  so  no  overestimate  occurs 

dA2up=dA2up+ (c2-xl ) *y ; 

121-121+sqrt ( (c2-xl) "2+ (y-yl) *2) ; 
else 

dA2up=dA2up+ (x2-xl) *y ; 

121=121+sqrt ( (x2-xl) ~2+ (y2-yl) ~2) ;  '/.Otherwise  calculate  area  and  add  to  sum 

end 


u=u+l; 

end 


u=u-l ; 
xcupper(u) ; 
dA2up 
121 


*/, Reset  counter  so  that  first  index  starts  where  we  left  off 
'/This  is  the  value  of  the  x  coordinate  where  we  left  off  above 
'/.Display  result  if  so  desired 
'/.Display  result  if  so  desired 


'/.Calculate  main  torque  box  lower  area  (dA21ow)  and  lower  skin  length  (122) 
dA21ow=0.0; 

122=0.0; 


while  xclower (j) <=c2 
xl=xclower (j ) ; 
x2=xclower (j  +  1) ; 
yl=abs (zclower ( j ) ) ; 
y2=abs (zclower ( j+1) ) ; 
x=xl+(x2-xl)/2; 

[y] =interp (xl , x2 , x , yl , y2) ; 


if  x2>c2 

dA21ow=dA21ow+(c2-xl)*abs (y) ; 
122=122+sqrt((c2-xl)“2+(y-yl)“2) ; 
else 

dA21ow=dA21ow+(x2-xl)*abs (y) ; 
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122=122+sqrt ( (x2-xl) ~2+ (y2-yl) ~2) ; 


end 

end 

'/,j=j-l;  '/.Reset  counter  so  that  first  index  starts  where  we  left  off 

'/.The  total  area  is  the  some  of  the  upper  and  lower  areas 

A2=dA2up+dA21ow 

’/********** ******************************************************** *********** 

’/.Calculate  trailing  edge  torque  box  upper  area  (A3)  and  upper  skin  length  (131) 

dA3up=0 . 0 ; 

131=0.0; 

while  u<=nupper-l 
xl=xcupper (u) ; 
x2=xcupper (u+1) ; 
yl=zcupper (u) ; 
y2=zcupper (u+1) ; 
x=xl+(x2-xl)/2; 

[y] =interp(xl ,x2,x,yl ,y2) ; 
dA3up=dA3up+ (x2-xl) *y ; 

131=131+sqrt ( (x2-xl) “2+ (y2-yl) *2) ; 
u=u+l ; 

end 

'/.Calculate  trailing  edge  torque  box  lower  area  (A3)  and  lower  skin  length  (132) 

dA31ow=0 .0; 

132=0.0; 

while  j<=nlower-l 
xl=xclower ( j) ; 
x2=xclower ( j+1) ; 
yl=abs (zclower (j) ) ; 
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y2=abs (zclower (j+1) ) ; 
x=xl+(x2-xl)/2; 

[y] =interp (xl , x2 , x , y 1 , y2) ; 

'/.In  case  of  reflexed  lower  surface:  If  lower  surface  below  chord  y<0,  and 
'/.incremental  area  is  added  (-*-=+).  If  lower  surface  above  chord  y>0,  and 
'/.incremental  area  is  subtracted  (-*+=-) 
dA31ow=dA31ow-(x2-xl)*y; 

132=132+sqrt ( (x2-xl) ~2+ (y2-yl) ~2) ; 

end 

dA31ow;  '/.Display  result  if  so  desired 
132  '/.Display  result  if  so  desired 

A3=dA3up+dA31ow 

y** ********* *************** 

'/.Find  spar  heights  (hi  and 
m=l ; 
n=l ; 

'/.Start  with  hi: 
while  xcupper (m) <=cl 

if  xcupper (m) <=cl  &  xcupper (m+1) >=cl 
xl=xcupper (m) ; 
x2=xcupper (m+1) ; 
yl=zcupper (m) ; 
y2=zcupper (m+1) ; 

[yll] =interp(xl ,x2, cl ,yl ,y2) ; 

end 

m=m+l ; 

end 

while  xclower (n) <=cl  '/.Find  main  spar  location  (=cl)  on  lower  xc 

if  xclower (n)  <=cl  &  xclower (n+1)  >=cl  '/.Find  if  around  main  spar  location 
xl=xclower (n) ; 


’/.Find  main  spar  location  (=cl)  on  upper  xc 
'/.Find  if  around  main  spar  location 


'/.Interpolate  for  upper  spar  height 


**************************************************** 

h2) 

'/.Upper  xc  loop  variable  initialized 
*/, Lower  xc  loop  variable  initialized 
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x2=xclower (n+1) ; 
yl=zclower (n) ; 
y2=zclower (n+1) ; 

[yl2]=interp(xl ,x2,cl ,yl,y2) ; 

end 
n=n+l ; 

end 

yii 

yi2 

hl=yll+abs(yl2) 

'/.Continue  with  h2: 

p=i; 

q=i; 

while  xcupper (p) <=c2 

if  xcupper (p)<=c2  &  xcupper (p+1) >=c2 
xl=xcupper (p) ; 
x2=xcupper (p+1) ; 
yl=zcupper (p) ; 
y2=zcupper (p+1) ; 
[y21]=interp(xl,x2,c2,yl ,y2) ; 

end 

p=p+l; 

end 

while  xclower (q) <=c2 

if  xclower (q) <=c2  &  xclower (q+1) >=c2 
xl=xclower (q) ; 
x2=xclower (q+1) ; 
yl=zclower (q) ; 
y2=zclower (q+1) ; 

[y22] =interp(xl ,x2 ,c2,yl ,y2) ; 

end 

q=q+l; 


'/, Interpolate  for  lower  spar  height 


'/, Calculate  main  spar  height 

'/.Upper  xc  loop  variable  initialized 
*/, Lower  xc  loop  variable  initialized 
’/.Find  rear  spar  location  (=c2)  on  upper  xc 
'/.Find  if  around  rear  spar  location 


*/,  Interpol  ate  for  upper  rear  spar  height 


'/.Find  rear  spar  location  (=c2)  on  lower  xc 
'/.Find  if  around  rear  spar  location 


'/.Interpolate  for  lower  rear  spar  height 
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end 


y21  '/.Display  result  if  so  desired 

y22  '/.Display  result  if  so  desired 

h2=y21+abs  (y22)  ’/.Calculate  rear  spar  height 

y#**********t*t************************************************t***t*************** 
return 
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function  [xcupper,zcupper,xclower,zclower]=NLF0215(xcupper  .zcupper  .xcloner  .zclower) 

7, Written  by  Capt  Peter  Cseke,  Jr.  -  Fall  1999  -  AFIT/ENY 
7, Defines  the  geometric  locus  of  the  airfoil  profile 
7.1nputs:  none 

7.0utputs :  none 

7.The  nondimensional  coordinates  of  the  NLF0215  high  altitude,  high  endurance  airfoil 


xcupper= [0 . 0  0.0024  0.00909  0.02004  0.03527  0.05469  0.07816 
0.10546  0.13635  0.17050. . . 

0.20758  0.24720  0.28894  0.33237  0.37702  0.42253  0.46864  0.51524  0.56247  0.61010. 

0.65752  0.70408  0.74914  0.79206  0.83222  0.86902  0.90193  0.93044  0.95409  0.97285. 

0.98710  0.99658  1.00000]; 

zcupper= [0 . 0  0.00917  0.01947  0.03027  0.04120  0.05201  0.06250 
0.07247  0.08175  0.09019... 

0.09761  0.10389  0.10887  0.11240  0.11428  0.11427  0.11219  0.10784  0.10147  0.09373. 

0.08513  0.07603  0.06673  0.05746  0.04844  0.03983  0.03175  0.02428  0.01737  0.01082. 

0.00507  0.00126  0.00000]; 

xclower= [0.00000  0.00245  0.01099  0.02592  0.04653  0.07242  0.10324 
0.13854  0.17788. . . 

0.22073  0.26654  0.31473  0.36468  0.41576  0.46731  0.51867  0.56920  0.61825  0.66662. 

0.71614  0.76645  0.81565  0.86198  0.90359  0.93862  0.96588  0.98504  0.99630  1.0]; 

zclower= [-0.00006  -0.00704  -0.01211  -0.01656  -0.02052  -0.02399 
-0.02699  -0.02954. .  . 

-0.03166  -0.03334  -0.03456  -0.03531  -0.03554  -0.03519  -0.03415  -0.03225  -0.02925 
-0.02441  -0.01663  -0.00705  0.00167  0.00804  0.01155  0.01198  0.00990  0.00655... 
0.00323  0.00086  0.0]  ; 


return 
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function  [y] =interp(xl ,x2 ,x,yl ,y2) ; 

'/.Written  by  Capt  Peter  Cseke,  Jr.  -  Fall  1999  -  AFIT/ENY 
'/.Linear  interpolation  subroutine 

*/, Inputs:  Maximum  and  minimum  independent  variables  x2  and  xl 

'/,  Maximum  and  minimum  dependent  variables  y2  and  yl 

*/,  Independent  variable  (x)  for  which  the  dependent  variable  (y)  is  sought 

'/.Outputs :  Value  y 

y=yl+(x-xl)*(y2-yl)/(x2“xl) ; 

return 
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'/, program  plotprof  lie  .m 


'/.Written  by  Capt  Peter  Cseke,  Jr.  -  Fall  1999  -  AFIT/ENY 
'/.Calculates  and  plots  the  geometric  properties  of  the  airfoil  profile 
'/.Uses :  NLF0215.m  subroutine  for  airfoil  input  data 

'/.Outputs:  plots  dimensional  and  non-dimensional  profile 

xcupper=[]  ;zcupper=[]  ;xclower=[]  ;zclower=[]  ; 

'/.Call  the  profile  function  for  non-dimensional  data 

[xcupper ,zcupper ,xclower .zclower] =NLF0215 (xcupper ,zcupper ,xclower .zclower) ; 

*/, These  are  some  of  the  design  variables: 

b=21 . 33; 

lam=0 .  40 ; 

croot=5.08; 

ctip=lam*croot ; 

tipoff set=(croot-ctip)/2; 

'/.Calculate  the  root  and  tip  chord  length  from  non-dimensional  data: 

'/,x=(x/ chord)* chord)  ;  z=(z/chord)  *chord)  ; 
xrootup=xcupper*croot ; 
zrootup=zcupper*croot ; 
xrootlo=xclower*croot ; 
zrootlo=zclower*croot ; 

'/.In  order  to  plot  the  tip  on  the  same  plot,  we  need  to  offset  it  according  to  taper 

xtipup=xcupper*ctip+tipof f set ; 

ztipup=zcupper*ctip ; 

xtiplo=xclower*ctip+tipof f set ; 

ztiplo=zclower*ctip ; 

'/.Plot  the  non-dimensional  profile  picture  from  airfoil  data 
close  all;  figure(l)  plot (xcupper .zcupper ,xclower ,zclower , Jb- ’ ) ; 
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title ( [’NLF(1)-0215F  Airfoil  -  Non-Dimensional  Coordinates’]); 
xlabel ( *  Station  (x/c) 7 ) ; 
ylabel( ’Ordinate  (z/c)’); 

axis ( [0  1  -0.5  0.5]); 
grid  on; 

’/.Plot  the  airplane-specific  root  and  tip  profiles: 
figure(2)  plot (xrootup,zrootup,xrootlo ,zrootlo , ’b-’) ; 
hold  on; 

plot (xtipup.ztipup, ’m- ’ ,xtiplo ,ztiplo , ’m- ’ ) ; 
hold  off; 

title ( [’NLF(1)-0215F  Airfoil  Root  and  Tip  Profiles’]); 
xlabel ( ’Station  (ft)’); 
ylabel ( ’Ordinate  (ft)’); 

axis([0  croot  -croot/2  croot/2]); 
grid  on; 
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Figure  A.l  The  NLF(1)-0215F  High-Endurance  Airfoil  Profile  at  Wing  Root  and  Tip 
(No  Washout) 
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Appendix  B .  Isotropic  Beam  Torsion 


Codes 


'/.program  shearflow3cell 

'/.Calculate  the  moments  required  and  shearflows  generated  by  achieving  desired 
'/.twist  per  unit  length  of  a  single-,  and  three-cell,  closed-section  torque  box  : 

*/. Written  by  Capt  Peter  Cseke,  Jr.  -  Fall  1999  -  AFIT/ENY 

'/, Calculates  the  shear  flow  in  three-cell,  closed  section  box  beam  due  to  imposed  tip 
*/, twist  angle.  It  also  calculates  the  necessary  concentrated  tip  moment. 

’/.Uses  :  area.m 

'/.Outputs:  shear  flows,  moment 

close  all; 
clear  all; 


*/, Input  Profile  parameters 

L=21 . 33 ; 

sparl=0. 25; 

spar2=0 .70 ; 

chord=5 .08; 


'/.Spanwise  lenght  of  torque  box  (ft) 
'/.Location  of  Main  Spar  ('/,  chord) 
'/.Location  of  Rear  Spar  ('/.  chord) 
'/.Chord  Length  at  station  (ft) 


plotprofile (L, chord)  ;  '/.Plot  the  picture  of  the  2-D  profile  (if  so  desired) 


tl=0.048; 
t2=tl ; 
t3=tl ; 
tsl=0 . 072 ; 
ts2=0 . 049 ; 


'/. Thickness 
’/.Thickness 
'/.Thickness 
'/.Thickness 
'/.Thickness 


(in)  of  nose  cone  skin  (Area  1) 
(in)  of  torque  box  skin  (Area  2) 

(in)  of  rear  surface  skin  (Area  3) 

(in)  of  main  spar 

(in)  of  rear  spar 
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c=chord*12 ; 
Lt=L*12; 


'/.Convert  chord  length  to  inches 
'/.Convert  spanwise  length  to  inches 


'/.Call  the  function  that  calculates  the  profile  geometric  parameters.  Send  spar 
'/.locations  as  non-dimensional  (x/c),  send  chord  length  (c)  in  inches. 

[A1 , 11 , A2 ,121 ,122 , A3 , 131 , 132 , hi , h2] = area (spar 1 , spar2 , c) ; 

A2=188 .4158; 

+*++****+++*+***********♦******+****************************************************** 
'/.Select  the  material  properties  of  the  beam: 

E=10 . 4*10~6;  '/.Young's  Modulus  for  Aluminum  2024-T3  [psi] 

G=3.86*10~6;  '/.Shear  Modulus  for  Aluminum  2024-T3  [psi] 

rhosp=0.101;  */, Density  for  Aluminum  2024-T3  [lb/in~3] 

Gsparl=G; 

'/.Designate  the  Shear  Modulus  of  the  spar  as  the  reference  shear  modulus 
G0=Gsparl ; 

'/.The  rear  spar  is  made  of  identical  material  as  main  spar: 

Gspar2=Gsparl ; 
rhosk=rhosp; 

Gskin=Gsparl ; 

'/(*  ************************************************************************* 

Gstar=Gskin/ GO ;  '/.Define  the  Weighted  Shear  Modulus 


'/************************************************************************** 
'/.Calculate  and  plot  moment  required  for  generating  the  given  angle  of  twist 
'/.for  single-cell ,  closed  section  torquebox. 

close  all; 

d=(l/2) *c* (spar2-sparl)  '/.Half  distance  between  spars 

thetamin=2;  '/.The  given  degrees  of  twist 
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thetamax=10; 

increment=2; 

j=0; 

'/.Generate  plot  of  required  moment  versus  angle  of  twist  for  single  cell  beam 
for  i=thetamin:  increment :  thetamax; 
theta=i; 

thetarad=theta*pi/ 180 ;  '/.The  twist  in  radians 

thetarad_unit=thetarad/Lt ;  '/.The  twist  per  unit  length  in  radians/in 

'/.Calculate  shear  flow  due  to  moment 

q(j)=thetarad_unit* (2*A2*G) * (1/ (121/t2+h2/ts2+122/t2+hl/tsl) ) 

'/.Calculate  the  moment  [in-lb]  required  for  twist: 

Mt=2*A2*q( j) ; 

M( j)=Mt/12 
angle ( j ) =theta ; 

end 

'/.Calculate  the  Moment  and  shear-flow  for  simplified  torquebox: 

T=sqrt  (l+(  (h2-hl)/(4*d))  "2)  ;  ’/.The  path  length  correction  constant 
A_s=2*d*  (0 . 5*hl+0 . 5*h2)  ;  '/.The  simplified  torque-box  area 

q_s=thetarad_unit* (2*A2*G) * (1/ (4*T*d/t2+h2/ts2+122/t2+hl/tsl) ) ; 

Mt_s=2*A_s*q_s ;  '/.The  Moment  (in-lb)  required  for  twist 

M_s=Mt_s/12;  '/.The  Moment  (ft-lb)  required  for  twist 

y>+  +  jjc  +  ****************3|c**3jc*  +  +  ****:************s»c***sjc*3|t*  +  ************************** 

'/.Calculate  cind  plot  moment  required  for  generating  the  given  angle  of  twist 
'/.for  three-cell,  closed  section  torquebox. 

'/.Calculate  the  constants  obtained  by  solving  the  Thetal  continuity  relations 
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/for  ql  (shear  flow  in  nose  torque  box)  in  terms  of  q2  and  q3,  and  substituting 
'/.back  into  the  Theta2=Theta3  continuity  relation,  solving  for  q3  (shear  flow  in 
/.trailing  edge  torque  box)  in  terms  of  q2: 

Cl=121/(Gstar*t2)+(h2/ts2)+(122/(Gstar*t2))+(hl/tsl)+(A2*h2/(A3*ts2)) ; 
C2=(ll/(Gstar*tl)+(hl/tsl)+(Al*hl/ (A2*tsl))) ; 

C3= (A2/A3) * (131/ (Gstar*t3) +132/ (Gstar*t3) +h2/ts2) ; 

C4= (A1/A2) * (121/ (Gstar*t2) +h2/ts2+122/ (Gstar*t2) +hl/tsl) ; 

/Substitute  q3  into  the  ql  equation,  and  solve  for  ql  in  terms  of  q2  only. 
/.Express  in  terms  of  constants: 

S3=(C1*C2- (C4*hl/tsl) -(hl/tsl) *2)/ (C2*C3+h2/ts2-(hl*h2/ (tsl*ts2) ) ) ; 

S1=(C4+ (hl/tsl) -(h2/ts2)*S3) / (11/ (Gstar*tl)+hl/tsl+(Al*hl/ (A2*tsl) ) ) ; 

**************************************************************************** 

'/, Calculate  the  effective  torsional  Stiffness  for  each  subsection 
GJeff 1=(4*A1~2)/ ( (11/ (Gskin*tl) )+(hl/ (GO*tsl) ) ) ; 

GJeff 2= (4* A2' 2) / ( (121/ (Gskin*t2) ) + (122/ (Gskin»t2) ) + (hi/ (G0*tsl) )+ (h2/ (G0*ts2) ) ) ; 
GJeff 3=(4*A3~2) / ( (131/ (Gskin*t3) ) + (132/ (Gskin*t3) ) + (h2/(G0*ts2) ) ) ; 

GJeff=GJeff l+GJeff2+GJeff3 

thetamin=2;  /, Degrees  of  twist  required  at  tip  (end  of  torquebox) 

j=0;  for  i=thetamin: increment : thetamax; 
theta=i ; 

j=j+i; 

angle (j)= theta; 

thetarad=theta*  (pi/180)  ;  ’/.Convert  degrees  to  radians 

thetaunit=thetarad/Lt ;  ’/, Radians  of  twist  per  unit  length  (rad/in) 

’/, The  applied  moment  required  to  achieve  the  required  twist  per  unit  length: 
Mt  (  j )  =  (GJeff  )  *thetaunit ;  */,in-lb 

’/.Then  the  shear  flow  in  the  second  subsection 
q2(j)=Mt(j)/(2*(Al*Sl+A2+A3*S3)) 
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'/Express  ql  and  q3  in  terms  of  the  constants  multiplying  q2: 
ql  (  j )  =Sl*q2  (  j )  '/,=ql=Mt/  (2*  (A1+A2/S1+A3*S3/S1) ) 

q3 ( j ) =S3*q2 ( j ) 

Mtftlb  ( j ) =Mt ( j ) /12  7,f  t-lb 

'/.The  total  moment  can  also  be  calculated  from: 

Mtsum=2*  (Al*ql+A2*q2+A3*q3)  ;  */, in-lb 

'/.This  should  yield  identical  answer  to  Mtftlb: 

Mt2=Mtsum/12 ;  '/.ft-lb 

'/.The  in-plane  shear  stress  (psi)  in  the  skin  of  the  nose  torque  box: 
sigmalskin(j)=ql (j)/tl; 

'/.The  in-plane  shear  stress  (psi)  in  the  skin  of  the  main  torque  box: 
sigma2skin ( j ) =q2 ( j ) /t 2 ; 

'/.The  in-plane  shear  stress  (psi)  in  the  skin  of  the  trailing  edge  box: 
sigma3skin(j)=q3(j)/t3; 

end 

'/.Plot  the  single-cell  and  three-cell  torque  box  required  moments  to  generate 
*/, given  angle  of  tip  twists 
figure(l) ; 

plot (angle ,M, *b~ * , angle , Mtftlb , *rs * ) ; 
hold  on; 

title ( [J Required  Moments  for  Single-  and  Three-Cell  Beam 
Torsion*] ) ; 

ylabeK* Moment  (ft-lb) ;); 
xlabel ( *Twist  Angle  (deg)*); 

legend( *  Single-Cell  Beam* , *  Three-Cell  Beam* ,4) ; 
plot (angle, M, ’b-* , angle, Mtftlb, *r-*) ; 
grid  on;  hold  off; 

'/.Plot  the  single-cell  and  three-cell  torque  box  shear  flows  generated  by 
'/.given  angle  of  tip  twists 
figure (2) ; 
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plot ( angle, q, 7b~7 , angle, ql, 7rs7 , angle, q2, ’go7 , angle, q3, 7md7) ;  hold 
on;  title ( [7 Shear  Flows  for  Single-  and  Three-Cell  Beam 
Torsion7]);  ylabel ( 7 Shear  Flow  (lb/in)7);  xlabel( 7 Twist  Angle 
(deg)7);  legend ( 7 Single-Cell  Beam7 , 7 Three-Cell  Beam,  Nose 
Section7 , . . . 

7Three-Cell  Beam,  Mid  Section7 , 7 Three-Cell  Beam,  Tail  Section7, 2); 
plot (angle ,q, 7b-7 , angle, ql, 7r-7 , angle, q2, 7g-7 , angle, q3, 7m-7) ;  grid 
on;  hold  off ; 
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Appendix  C.  Anisotropic  Composite 


Beam  Torsion  Codes 

'/.Program  Compl_aniso_compz  .m 

*/, Written  by  Capt  Peter  Cseke,  Jr.  -  Fall  1999  -  AFIT/ENY 
'/.Driver  program  to  Matlab  functions  Layup. m 
%  Anisotorsion_compz  .m 

7, Defines  trapezoid  single-cell  box  beam  geometric  and  material  properties,  applied 
'/.cross  section  loads  and  moments.  Calculates  tip  twist  angle  of  box  beam. 

7. 

'/.Uses  :  Layup. m  to  calculate  engineering  properties  of  composite  lamina. 

'/,  Anisotorsion_compz .m  to  calculate  twist  angle  of  box  beam. 

'/.Outputs:  Tip  twist  angle  of  box  beam. 

7. 

'/.User  enters : 

'/.orientation:  [0  45  -45  90  -90]  etc,  as  a  row  vector; 

'/.times  :  The  number  of  times  ply  is  repeated; 

'/.symmetry  :  0  for  no; 

7.  1  for  yes 

'/.Example  :  [0  45  -45  90]  4s 

7,  orient  -  [0  45  -45  90] 

7,  times  =  4 

'/,  sym  =  1 

'/.Note  :  all  properties  are  vectors  with  properties  per  ply. 

7, Limits  :  Assumes  weighted  center  of  laminate  is  at  h/2. 

7#* ******************************************************************************* 
close  all; 
clear  all; 

'/.Define  beam  geometric  parameters 

cr=60.96;  '/.Root  chord  length  [in] 
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cl=0.25;c2=0.70; 
hl=8.4334; 
h2=5 . 3035 ; 
121=27.7732 ; 
122=29.2992; 
L=255 . 96; 


’/.Main  spar  and  Rear  spar  location  in  chord  percent 
'/.Height  if  Main  Spar  [in] 

'/, Height  of  Rear  Spar  [in] 

'/.Length  of  top  surface  skin  between  Main  and  Rear  Spars  [in] 
'/.Length  of  bottom  surface  skin  between  Main  and  Rear  Spars  [in] 
'/.Spanwise  length  of  torque  box  [in] 


'/.Define  the  Loads  and  Moments  acting: 


P=0; 

Vx=0 ; 

Vy=0; 

Mx=0; 

My=0; 

M_ftlb=4404; 

M=M_ftlb*12; 


'/, Spanwise  extensional  load  [lb] 

'/, Vertical  shear  load  (lift)  [lb] 

'/, Horizontal  shear  load  (drag)  [lb] 

'/.Applied  moment  about  x  axis  (bending)  [in-lb] 
'/.Applied  moment  about  y  axis  (bending)  [in-lb] 
*/, Applied  moment  about  z  axis  (torque)  [ft-lb] 

7.  [in-lb] 


y* **************************************************** 
'/.Carbon/Epoxy  (AS4/3501-6) 


El=20.6*10~6; 
E2=l.  50+1CT6; 
G12=l .04*10~6 ; 
vl2=0 . 27 ; 
v21=0 . 02; 


7, Youngs  Modulus  (fiber  direction) 
'/.Young's  Modulus  (matrix  direction) 
'/.Shear  Modulus  (in-plane) 

'/.Poisson's  Ratio 
'/.Poisson's  Ratio 


y#* ********************************************************************************* 


theta=45; 
times=l ; 
sym=l ; 

thick=0 . 005 ; 
t=  [thick] ; 


'/.Center  composite  fiber  rotation  angle  (deg) 

'/.Number  of  time  lamina  is  wound 
'/.Symmetry  of  lamina  0=no,l=yes 

'/.Thickness  per  layer  of  angle  lamina  (given  by  material  props.) 
’/.Thickness  vector  per  layer  of  lamina 


y* ********************************************************************************** 
'/.Define  layer  of  angled  lamina  (only  one  layer  of  variable)  orientation  angles 
anglemin=10; 
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angle incr=5; 
anglemax=15; 

count=0;  7, Reset  the  loop-count  variable  (#  of  angles) 

j=i; 

for  j=anglemin: angle incr: anglemax  7, Loop  through  the  angle  range 

count=count+l ; 

angle  (count)  =j  ;  7.Form  the  angle  vector  for  the  plot 
’/,******  ******  ***************************************  ******************************** 
7.Define  one  layer  of  angled  lamina  orientation  angles 
orientcom= [45  45  -45  -45  0  0  90  90]  ; 

y*********************************************************************************** 

7, Redefine  thicknesses  according  to  number  of  plies 

tl=times*thick*length(orient)  *  (sym+1)  ;  7iTotal  thickness  of  laminate  at  main  spar 

t2=tl;  7.Total  thickness  of  laminate  at  rear  spar 

ts=tl;  7.Total  thickness  of  laminate  at  skin 

y* ********************************************************************************** 
7.Form  the  material  properties  for  each  angle  ply  (number  of  angle  ply  from  orient) 
Elx=ones (size (orientcom) ) *E1 ; 

E2x=ones (size (orientcom) ) *E2 ; 

G12x=ones (size (orientcom) ) *G12 ; 
vl2x=ones (size (orientcom) ) *vl2 ; 
tx=ones (size (orientcom) ) *t ; 

7.Calculate  engineering  properties  of  the  composite  laminate  in  the  structural  axes : 

[Ex ,Ey , Gxy , a , vxy , vyx , nsx ,nxs , nys , nsy] =Layup (Elx , E2x , vl2x , G12x , tx , orientcom , times , sym) ; 

7.Calculate  the  tip  twist  angle  for  the  composite  construction  torquebox 

[S,f ideg]=Anisotorsion_compz(cr,cl , c2,hl ,h2 ,121 ,122,tl ,t2,ts ,L,P, Vy,Vx,Mx,My ,M, . . . 

Ex,Ey ,Gxy ,vxy , vyx,nsx,nxs ,nys ,nsy) ; 
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twist (count) =f ideg ; 


end 

twist 

'/.Declare  plot  style  lines  and  markers 
linestyle=['b-  ' ; ' g-  '  ; 'r-  ' ; 'm-  ';'c-  9  ; *k-  1 ; 'y-  '] 
markstyle=['b+  ';'go  ';'rh  ' ; 'md  *;'cp  ' ; 'k*  ' ; 'yv  '] 

figure (1)  plot (angle, twist, markstyle (1, :)) ;  hold  on; 
plot (angle, twist, linestyle ( 1, :)) ;  hold  on; 

title ( ['Twist  Angles  Due  To  Uniform  Moment']); 

ylabel ( ' Twist  Angle (deg) ' ) ; 
xlabel( 'Fiber  Angle  (deg)'); 
grid  on; 
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'/.Function  Anisotorsion.compz 

7, Written  by  Capt  Peter  Cseke,  Jr.  -  Fall  1999  -  AFIT/ENY 

'/.Calculates  the  tip  twist  angle  of  the  single-cell  trapezoid  box-beam  using 
'/.Libove’s  solution 

'/.Uses :  Box  beam  geometric  and  material  parameters  are  input  from  driver  program 

'/.Outputs:  Tip  twist  angle  in  degrees. 

*/, Comment:  Also  works  for  isotropic  torsion,  due  to  fact  that  Beta2=0  for 

'/,  isotropic  materials. 

function 

[S,Fideg]=Anisotorsion_compz(cr,cl,c2,hl,h2,121,122,tl,t2,ts,L,P,Vy,Vx,Mx,My,M, . . . 

Ex , Ey , Gxy , vxy , vyx , nsx , nxs , ny s , nsy ) ; 

y#%+^+t***********************************************************+********** 

'/.Define  the  half-width  (d)  ,  and  top  surface  length  correction  factor  (T) 
d=cr* (c2-cl) /2 ; 

T=sqrt(l+((h2-hl)/(4*d))~2) ; 

Area=2*d* (hl/2+h2/2) ; 

'/.Calculate  distance  from  origin  to  top  and  bottom  surfaces  (approximation) 
m=(hl+h2)/4; 

'/.The  Compliance  Matrix  in  the  x-y  structural  axes: 

Sll^l/Ex;  S12=-vyx/Ey;  S14=nsx/Gxy; 

S21=-vxy/Ex;  S22=l/Ey;  S24=nsy/Gxy; 

S41=nxs/Ex;  S42=nys/Ey;  S44=l/Gxy; 

S=[S11  S12  S14; S21  S22  S24;S41  S42  S44] ; 

'/.Define  elastic  constants  for  constitutive  relations: 
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alf al=Sll ; alf a2=S14 ; alf a4=S44 ; 


•/.Define  elastic  constants  for  force-strain  relations: 
Bl=l/alfal;B2=-alfa2/alfal;B4=alfa4-(alfa2~2/alfal); 

7, Define  the  elements  of  the  B  matrix: 

bll=4*Bl*ts*d*T+Bl*t2*h2+Bl*tl*hl;  bl2=0;  bl3=-Bl*d* (t2*h2-tl*hl) ; 

b21=0 ; 

b22=(-l/12)  *B1* (t2*h2~3+4*ts*d*T*h2~2+4*ts*d*T*hl*h2+tl*hl~3+4*ts*d*T*hl~2)  ; 
b23=0; 

b31=Bl*d* (t2*h2-tl*hl) ; 
b32=0; 

b33=(-l/3)*Bl*d~2>K4*ts*d*T+3*t2*h2+3*tl*hl) ; 

B=[bll  bl2  bl3 ;b21  b22  b23;b31  b32  b33] ; 

A=inv(B) ; 

all=A(l,l) ;al2=A(l,2) ;a!3=A(l,3) ; 
a21=A(2,l) ;a22=A(2,2) ;a23=A(2,3) ; 
a31=A(3, 1) ; a32=A(3, 2) ; a33=A(3,3) ; 

a4=8*m*Bl*d~2*T''2*ts+4*d~2*Bl*hl*ts*T+d*Bl*hl*t2*h2+  (1/2)  *d*Bl*tl*hl~2+ .  .  . 
4*d~3*Bl~2*T~2*t2~2*h2+2*d~2*Bl*T*ts*h2+d~2*Bl~2*T*t2~2*h2~2+ .  .  . 
2*m*Bl*d*T*t2*h2; 

a5=(l/12)  *Bl*d*  (-tl*hl''3+t2*h2~3+6*ts*T*d*h2~2+6*ts*T*d*h2*hl+8*m*ts*T''2*d*h2+ . 
16*m*ts*T~2*d*hl)  ; 

a6=(l/2)*Bl*d~2*(4*m*t2*T*h2+2*hl*t2*h2-tl*hl~2+t2*h2~2) ; 

/.Define  the  elements  of  the  C  matrix: 
cll=B2*h2-B2*hl ;  c21=B2*d*T* (h2+hl) ;  c31=B2*d* (h2+hl) ; 

cl2=4*Bl*ts*d~2*T''2*B2-4*B2*Bl~2*d~2*T~2*t2~2*h2-2*B2*Bl*d*T*ts*h2 .  .  . 
-B2*Bl~2*d*T*t2~2*h2~2+2*B2*Bl*d*T*t2*h2+4*Bl*B2*hl*ts*d*T . . . 
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+B2*Bl*hl*t2*h2+B2*Bl*hl~2*tl/2; 

c22=(-l/12)*B2*Bl*(-tl*hl~3+6*h2~2*t2*d*T+6*h2*hl*t2*d*T+24*ts*h2*d~2*T''2- . . . 
2*d*T*Bl*t2~2*h2''3) ; 

c32=(-l/6) *B2*Bl*d* (8*ts*d~2*T~2+6*hl*h2*t2+12*ts*T*d*h2+24*ts*hl*d*T . . . 
+3nl*hl~2+24*T~2*d~2*Bl*t2~2*h2+6*d*T*Bl*t2~2*h2~2) ; 

cl3=(l/12)*Bl*B2*(tl*hl~3+6*ts*d*T*h2~2+6*ts*d*T*hl*h2+t2*h2~3) ; 

c23=  (1/4)  *  (Bl*B2*ts*T~2*d''2*  (h2~2+hl~2+2*hl*h2) )  ; 

c33= (1/12) *B1*B2* (-tl*hl~3+6*ts*d*T*h2~3+6*ts*d*T*hl*h2+t2*h2~3) ; 

cl4=(“l/6)*Bl*B2*d*(-3nl*hl~2+6*hi*h2n2+8*ts*d~2*T~2+12*t2*d*T*h2-3*t2*h2~2) ; 
c24=  (1/12)  *Bl*B2*d*  (6*h2~2*t2*d*T+6*h2*hl*t2*d*T-h2~3*t2+tl*hl''3)  ; 
c34=(l/2) *Bl*B2*d~2* (2*hl*h2*t2“tl*hl~2+h2~2*t2) ; 


C=  [ell  cl2  cl3  cl4; c21  c22  c23  c24;c31  c32  c33  c34] ; 

*/, Define  the  elements  of  the  D  vector: 
dl=B4*(4*d*T*tl*t2+ts*tl*h2+ts*t2*hl)/(ts*t2*tl) ; 
d2=(-l/2)*B4*Bl*(8*d~2*T~2*ts*tl*Bl*t2~2*h2+4*d*T*ts~2*tl*h2. . . 

+2*d*T*ts*tl*Bl*t2~2*h2~2+4*t2~2*tl*d*T*h2+16*d~2*T~2*t2*ts*tl+. . . 
8*t2*ts~2*hl*d*T+2*ts~2*ts*hl*h2+t2*ts*tl*hl~2)/(t2*ts*tl)  ; 
d3=(l/12)  *B4*B1*  (“t2*hl~3+t2*h2''3+6*ts*d*T*h2~2+6*ts*d*T*hl*h2+16*t2*d'‘2*T~2*hl .  .  . 
+8*t2*d''2*T~2*h2)  /t2 ; 

d4=(-l/2) *B4*Bl*d* (-4*tl*t2*d*T*h2-ts*tl*h2~2+ts*tl*hl~2-2*ts*hl*t2*h2) / (ts*tl) ; 

‘/.Define  the  shear  flow  at  s=0: 

qO=(l/(2*Area) ) * (M-a22*a5*Vy+ (al3*a4-a33*a6) *Vx) ; 

’/, Define  the  derivatives  of  the  longitudinal  strain,  and  curvatures: 
e0p=al3*Vx; 

Kxp=a22*Vy ; 

Kyp=a33*Vx; 

*/, Define  the  values  for  the  Q’s: 
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Ql=cll*qO+cl2*eOp+cl3*Kxp+cl4*Kyp; 

Q2=c21*q0+c22*e0p+c23*Kxp+c24*Kyp; 

Q3=c31*q0+c32*e0p+c33*Kxp+c34*Kyp; 


eO=all* (P-Ql) -al2* (Mx+Q2) -al3* (My+Q3) ; 

Kx=a21* (P-Ql) -a22* (Mx+Q2) -a23* (My+Q3) ; 

Ky=a31* (P-Ql) -a32* (Mx+Q2) -a33* (My+Q3) ; 

7, The  rotation  in  radians  per  unit  length: 

dFdz= (1/ (2*Area) ) * (-e0*cll+Kx*c21+Ky*c31+q0*dl+e0p*d2+Kxp*d3+Kyp*d4) ; 

7, The  rotation  in  degrees  per  unit  length: 
df dz_deg=dFdz*180/pi ; 

7»The  rotation  in  degrees  for  entire  length  of  beam: 

Fideg=df dz_deg*L ;  return 
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'/.Function  Layup 


*/, Calculates  engineering  properties  of  laminate  consisting  of  uniform  lamina, 
’/.and  orientation  angles . 

'/, Input :  Composite  lamina  properties  in  material  1-2  axes 

'/.Return:  Laminate  engineering  properties  for  the  purposes  of  stiffness  (C)  and 
'/,  compliance  matrix  (S)  calculations. 

function 

[Ex , Ey , Gxy , a , vxy , vyx , nsx , nxs , nys , nsy] =Layup (El , E2 , vl2 , G12 , t , orient , times , sym) 

y( +*+**♦*************************************************=♦'********************** 

'/.  Layup  (El,  E2,  vl2,  G12,  t,  orient,  times,  sym) 

%  El:  lamina  Young’s  Modulus  (1  direction)  per  ply 

*/,  E2:  lamina  Young’s  Modulus  (2  direction)  per  ply 

'/,  vl2:  lamina  Poisson’s  Ratio  per  ply 

'/.  G12:  lamina  Shear  Modulus  per  ply 

%  t:  vector  of  lamina  thickness  per  ply 

*/,  orient :  vector  of  orientations  in  degrees 

%  times:  Multiples  of  orient 

%  sym:  Symmetric?  (0  =  NO,  anything  else  means  yes) 

X 

%  Example:  [0  45  -45  90] 4s 

%  orient  =  [0  45  -45  90] 

*/,  times  =  4 

*/.  sym  =  1 

% 

*/,  note:  all  properties  are  vectors  with  properties  per  ply. 

'/,  Assumes  midplane  =  h/2. 

y. 

'/,  Written  by:  Capt  Jim  Rogers 

’/,  Modified  by:  Capt  Peter  Cseke 

A=zeros (3,3) ;  B=A;  D=A ; 
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for  i=l:times 


Elex((i-l)*length(El)+l : i*length(El) )=E1 ; 

E2ex ( (i-1) *length(E2) +1 : i*length(E2) ) =E2 ; 
vl2ex((i-l)*length(vl2)+l : i*length(vl2) ) =vl2 ; 

G12ex( (i-1) *length(G12) +1 : i*length (G12))=G12; 
tex ( (i-1) *length(t) +1 : i*length (t) ) =t ; 

orientex( (i-1) *length(orient)+l : i*length (orient) ) =orient ; 
end  if  syra  ~=  0 

Elex(2*length(Elex) :-l : length (Elex)+l)=Elex; 
E2ex(2*length(E2ex) : -1 : length(E2ex) +l)=E2ex; 
v!2ex(2*length(vl2ex) : -1 : Iength(vl2ex)+l)=vl2ex; 

G12ex (2*length (G12ex) :  -1 : length (G12ex) +1) =G12ex ; 
tex(2*length(tex) : -1 : length (tex) +1) -tex; 

orientex (2*length (orientex) : -1 : length (orientex) +1) =or ientex ; 

end 

'/.The  number  of  laminas  (in  one  ply) 
n=length (orientex) ; 

tott=sum(tex) ;  h^zeros (n+1 , 1) ; 

h(l)=-tott/2;  for  i=2:n+l 

h(i) =sum(tex(l : i-1) ) -tott/2 ; 

end 

for  i=l:n 

'/.From  the  symmetry  of  the  compliance  matrix: 
v21=vl2ex(i)*E2ex(i)/Elex(i) ; 
qll=Elex(i)/(l-vl2ex(i) *v21) ; 

'/.Equation  (3.56) 

Q12= [  qll  v21*ql 1  0; 

v21*qll  qll*E2ex(i)/Elex(i)  0; 
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0 


0 


G12ex(i)]  ; 


'/.Equat ion  (3.59) 

mc=cos (orientex(i) *pi/180) ; 

ns=sin(orientex(i)*pi/180) ; 

T=[mc~2  ns~2  2*mc*ns;  ns~2  mc~2  -2*mc*ns;  -mc*ns  mc*ns  mc~2-ns~2] ; 
Q12(:,3)=2*Q12(:,3); 

'/.Equation  (3.66) 

Qxy=T\Q12*T; 

'/.Reset  the  actual  value  of  the  3rd  column  of  Q12 
Q 12 ( : , 3) =Q 12 ( : ,3) /2; 

'/.Divide  3rd  column  by  2 
Qxy(: ,3)=Qxy(: ,3)/2; 

*/, Calculate  laminate  stiffness  matrices  (Eq.  5.20) 

A=A+Qxy*(h(i+l)-h(i) ) ; 

B=B+Qxy*(h(i+l)~2-h(i)~2)/2; 

D=D+Qxy* (h(i+l) ~3-h(i) ”3) /3; 

end 

'/.Calculate  the  laminate  compliance  matrices  (Eq.  5.27) 

Bs--A\B ;  Cs=B/A;  Ds=D-B* (A\B) ; 

'/.The  laminate  extensional  compliance  matrix  (S)  : 
a=inv (A) -Bs* (Ds\Cs) ; 

'/.Calculate  laminate  engineering  properties  (barred  values)  referenced  to 
*/,x  and  y  axes  (including  Poisson’s  ratios,  and  shear-coupling  coefficients) 
Ex=l/tott/a(l , 1) ; 

Ey=l/tott/a(2,2) ; 

Gxy=l/tott/a(3,3) ; 
vxy=-a(2 ,l)/a(l,l) ; 
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vyx=-a(l ,2) /a (2 ,2) 
nsx=a(l ,3) /a(3 ,3) ; 
nxs=a(3, 1) /a(l , 1)  ; 
nys=a(3,2) /a(2 ,2) ; 
nsy-a(2,3)/a(3,3) ; 


return 


Appendix  D.  Anisotropic 


Composite-PZT  Beam  Torsion  Codes 

'/.Program  Compl_aniso_pzt  .m 

*/, Written  by  Capt  Peter  Cseke,  Jr.  -  Fall  1999  -  AFIT/ENY 

*/, Driver  program  to  Matlab  functions  Layup. m 

'/.  Anisotorsion_pzt  .m 

'/.Defines  trapezoid  single-cell  box  beam  geometric  and  material  properties,  applied 
'/.cross  section  loads  and  moments.  '/.Calculates  tip  twist  angle  of  box  beam. 

•/. 

'/.Uses  :  Layup. m  for  engineering  properties  of  lamina  and  composite-pzt  lamina. 

'/,  Anisotorsion_pzt  .m  to  calculate  twist  angle  of  box  beam. 

'/.Outputs:  Tip  twist  angle  of  box  beam. 

y. 

'/.User  enters : 

'/.orientation:  [0  45  -45  90  -90]  etc,  as  a  row  vector; 

'/.times  :  The  number  of  times  ply  is  repeated; 

'/.symmetry  :  0  for  no; 

'/.  1  for  yes 

'/.Example  :  [0  45  -45  90]  4s 

'/.  orient  =  [0  45  -45  90] 

*/,  times  =  4 

'/,  sym  =  1 

'/.User  enters:  Loading  condition  on  cross  section  (P , Vx,V,Mx,My ,M=Mz) 

'/,  Select  PZT  laminate  material  properties  from  available  choices 

'/.Note  :  all  properties  are  vectors  with  properties  per  ply. 

'/.Limits  :  Assumes  weighted  center  of  laminate  is  at  h/2. 

•yf****  *************  ******  ********************************************************* 

close  all; 
clear  all; 
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'/.Define  beam  geometric  parameters 


cr=60 . 96; 

cl=0 . 25 ; c2=0 . 70; 

hl=8.4334; 

h2=5 . 3035 ; 

121=27.7732; 

122=29.2992; 

Xt 1-0. 072; 

Xt 2=0. 049; 

'/,ts=0 .048; 

L=255 . 96; 


'/.Root  chord  length  [in] 

'/.Main  spar  and  Rear  spar  location  in  chord  percent 
'/.Height  if  Main  Spar  [in] 

’/.Height  of  Rear  Spar  [in] 

’/, Length  of  top  surface  skin  between  Main  and  Rear  Spars  [in] 
'/.Length  of  bottom  surface  skin  between  Main  and  Rear  Spars  [in] 
'/.Main  Spar  thickness  [in] 

’/.Rear  Spar  thickness  [in] 

’/.Skin  thickness  (tl=t2)  [in] 

’/.Spanwise  length  of  torque  box  [in] 


'/.Define  the  Loads  and  Moments  acting: 


P=0; 

Vx=0 ; 

Vy=0; 

Mx=0; 

My=0; 

M_ftlb=4404; 
M=M_f tlb*12 ; 


X Spanwise  extensional  load  [lb] 
'/.Vertical  shear  load  (lift)  [lb] 
'/.Horizontal  shear  load  (drag)  [lb] 

'/, Applied  moment  about  x  axis  (bending) 
’/.Applied  moment  about  y  axis  (bending) 
'/.Applied  moment  about  z  axis  (torque) 

X 


[in-lb] 

[in-lb] 

[ft-lb] 

[in-lb] 


***************************************************** 
'/.Carbon/Epoxy  (AS4/3501-6) 


El=20. 6*1CT6; 
E2=l .  50+1CT6; 
G12=l. 04*10*6; 
vl2=0. 27 ; 
v21=0.02; 
thick=0 .005 ; 


'/.Young's  Modulus  (fiber  direction) 

*/, Young's  Modulus  (matrix  direction) 

*/, Shear  Modulus  (in-plane) 

'/.Poisson's  Ratio 
'/.Poisson's  Ratio 

'/.Thickness  per  layer  of  angle  lamina  (given  by  material  props.) 


x **************************************** ************** 

*/, Piezoceramic  Lamina 

Elpzt=4. 6786*10*6 ;  '/.Young's  Modulus  (AFC)  (poling,  long,  or  fiber  direction  psi) 
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E2pzt=2.4173*l(T6;  */, Young’s  Modulus  (AFC)  (transverse  direction)  (psi) 

G12pzt=5 . 8015*10~5 ;  ’/.Shear  Modulus  (AFC)  (in-plane ,  psi) 

thickpzt=0. 006496063;  '/, Thickness  of  one  layer  of  PZT  lamina  (AFC) 

'/,Elpzt=5.4389*lCT6;  '/.Young’s  Modulus  (Eon  Barrett) 

'/.E2pzt=2. 0305  +  1CT6;  '/.Young’s  Modulus  (Ron  Barrett) 

*/,G12pzt=5 . 5114*10*'5 ;  '/.Shear’s  Modulus  (Ron  Barrett) 

*/,thickpzt=0 . 005 ;  */, Thickness  of  one  layer  of  PZT  lamina  (Ron  Barrett) 

'/,dl3=l  . 66*10~-10;  ’/.Dielectric  constant  (m/V) 

'/,Elpzt=9. 1374*10~6;  '/.Young’s  Modulus  (Zhou,  Liang  &  Rogers) 

'/,Elpzt=10. 153*10~7;  '/.Young’s  Modulus  (Crawley,  de  Luis) 

vl2pzt=0.30;  '/.Poisson’s  Ratio 

*/,v21pzt=0 . 30;  ’/.Poisson ’s  Ratio 

yi+*+*****t********************************************:*t************************:f** 

anglepzt=45;  '/.Orientation  of  PZT  patch 

theta=45;  '/.Center  composite  fiber  rotation  angle  (deg) 

times=l;  '/.Number  of  times  lamina  is  wound 

sym=0;  */, Symmetry  of  lamina  0=no,l=yes 

y+^+t******************************************************************************* 

'/.Define  layer  of  angled  lamina  (only  one  layer  of  variable)  orientation  angles 

anglemin=0 ; 

angleincr=2; 

anglemax=90; 

count=0;  '/.Reset  the  loop-count  variable  (#  of  angles) 

for  j  =anglemin:  angle incr :  anglemax  '/.Loop  through  the  angle  range 

count=count+l ; 
angle (count) = j ; 
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'/(**  *******************************************************************  ************ 

'/, Define  layer  of  angled  lamina-pzt  orientation  angles 
orientall=[0  j  anglepzt  -j  0]; 
orientcom=[0  j  theta  -j  0]; 

locatepzt=3;  ‘/.The  PZT  layer  location  in  the  laminate 

yt********************************************************************************* 
'/.Redefine  thicknesses  according  to  number  of  plies 

t=  [thick];  '/.Thickness  vector  per  layer  of  lamina 

tl=times*thick*length(orientcom)  *  (sym+1)  ;  */, Total  thickness  of  laminate  at  main  spar 
t2=tl;  '/.Total  thickness  of  laminate  at  rear  spar 

'/.Total  thickness  of  laminate  (composite  and  pzt)  at  skin 
ts=times*thick* (length(orientall) -1) * (sym+1) +thickpzt ; 

y,* ****************************************+***************************************** 

'/.Form  material  properties  for  each  angle  ply  (number  of  angle  ply  from  orientall) 
Elxp=ones (size (orientall) ) *E1 ; 

E2xp=ones (size (orientall) ) *E2 ; 

G12xp=ones (size (orientall) ) *G12 ; 
vl2xp=ones (size (orientall) ) *vl2 ; 
txp=ones (size (orientall) ) *t ; 

Elxp (locatepzt) =Elpzt ; 

E2xp (locatepzt) =E2pzt ; 

G12xp (locatepzt) =G12pzt ; 
vl2xp (locatepzt) =vl2pzt ; 
txp (locatepzt) =thickpzt ; 

*/t*  ************  +  +  *********+******************************************************♦** 
'/.Call  Layup. m  subroutine  to  calculate  composite-pzt  laminate  engineering  properties 
[Epx.Epy , Gpxy , ap , vpxy , vpyx.npsx ,npxs ,npys ,npsy]=Layup (Elxp, E2xp, vl2xp,G12xp, txp, . . . 

orientall , times , sym) ; 


D-4 


’/.The  Compliance  Matrix  (plane  stress)  for  the  Composite-PZT  Laminate  in  the  x-y 
•/.structural  axes : 

Spll=l/Epx;  Spl2=-vpyx/Epy;  Spl4=npsx/Gpxy ; 

Sp21=-vpxy/Epx;  Sp22=l/Epy;  Sp24=npsy/Gpxy ; 

Sp41=npxs/Epx;  Sp42=npys/Epy ;  Sp44=l/Gpxy ; 

Spzt=  [Spll  Spl2  Spl4; Sp21  SP22  Sp24;Sp41  Sp42  Sp44] ; 

%************+********************************************************************* 
'/.Calculate  laminate  engineering  properties  (composite  laminate  only) 

Elx=ones (size (orientcom) ) *E1 ; 

E2x=ones (size (orientcom) ) *E2 ; 

G12x=ones (size (orientcom) ) *G12 ; 
vl2x=ones (size (orientcom) ) *vl2 ; 
tx=ones (size (orientcom) ) *t ; 

'/.Call  Layup. m  subroutine  to  calculate  composite  only  laminate  engineering  properties 
[Ex , Ey , Gxy , a , vxy , vyx ,nsx , nxs , nys , nsy] =Layup (Elx , E2x , v!2x , G12x , tx , orientcom .times , sym) ; 

'/.The  Compliance  Matrix  (plane  stress)  for  the  Composite-PZT  Laminate  in  the  x-y 
'/.structural  axes : 

Sll=l/Ex ;  S12=-vyx/Ey ;  S14=nsx/Gxy; 

S21=-vxy/Ex;  S22=l/Ey;  S24=nsy/Gxy; 

S41=nxs/Ex;  S42=nys/Ey;  S44=l/Gxy; 

Scom= [Sll  S12  S14; S21  S22  S24;S41  S42  S44] ; 

y(^+++*************************************************+:*******************t********** 
'/.Call  anisotropic  torsion  solution  subroutine  to  calculate  the  amount  of  tip  torsion 
'/.for  the  single  cell  box  beam 

[f ideg]  =Anisotorsion_pzt (cr , cl , c2 ,hl ,h2 ,121 ,122 , tl ,t2 ,ts ,L ,P , Vy , Vx.Mx.My ,M, Spzt , Scorn)  ; 
twist ( count )=fideg; 
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end 


orientall 

twist 

'/.Declare  plot  style  lines  and  markers 


linestyle= [7b-  3 

Vg“  ’ 

;’r-  ’ 

;7m-  7 ;  7c-  7 

;7k-  7 ; 7y- 

»] 

markstyle= [7b+  3 

;>g0  3 

;’rh  1 

; 7md  7 ; 7cp  7 

; 7k*  7;7yv 

’] 

figure (1) ;  plot (angle, twist ,markstyle (2, : )) ; 
hold  on; 

plot (angle , twist ,linestyle (2 , : ) ) ; 
hold  on; 

title ( [7 Twist  Angles  Due  To  Uniform  Moment7]); 

ylabel( 7 Twist  Angle  (deg)’); 
xlabel ( 7 Fiber  Angle  (deg)7); 
legend  ( 7  AFC  \0); 
grid  on; 
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'/, Function  Anisotorsion_pzt 


'/.Written  by  Capt  Peter  Cseke,  Jr.  -  Fall  1999  -  AFIT/ENY 

‘/.Calculates  the  tip  twist  angle  of  the  single-cell  trapezoid  box-beam  with  PZT  actuator 
‘/.embedded  in  top  and  bottom  skins.  No  PZT  in  main  or  rear  spars. 

‘/.Uses:  Box  beam  geometric  and  material  parameters  are  inputed  from  driver  program 

‘/.Output:  Tip  twist  angle  in  degrees 

‘/.Also  works  for  isotropic  torsion,  because  Beta2=0  (B2  and  Bp2)  for  isotropic  materials. 

function  [Fideg] =Anisotorsion_pzt (cr , cl , c2 ,hl ,h2 , 121 , 122 ,tl ,t2 ,ts ,L,P , Vy , Vx ,Mx,My ,M, Sp,Sc) ; 

’^* ************************************************************************** 

‘/.Define  the  half -width,  and  top  surface  length  correction  factor 
d=cr* (c2-cl) /2; 

T=sqrt (1+ ( (h2-hl) / (4*d) ) ~2) ; 

Area=2*d* (hl/2+h2/2) ; 
m=(hl+h2)/4; 

‘/.Designate  the  compliance  Matrix  for  the  composite  laminate: 

Sll=Sc(l,l) ;S12=Sc(l,2) ;S14=Sc(l,3) ; 

S21=Sc(2,l) ;S22=Sc(2,2) ;S24=Sc(2,3) ; 

S41=Sc(3,l) ;S42=Sc(3,2) ; S44=Sc (3,3) ; 

'/, Designate  the  compliance  Matrix  for  the  composite-pzt  laminate: 

Spll=Sp(l ,  1)  ; Spl2=Sp(l , 2) ; Spl4=Sp(l ,3) ; 

Sp21=Sp(2 , 1) ; Sp22=Sp(2 ,2) ; Sp24=Sp(2 ,3) ; 

Sp41=Sp (3,1); Sp42=Sp (3 , 2) ; Sp44=Sp (3 , 3) ; 

'/.Define  elastic  constants  for  constitutive  relations  of  composite  laminate  only: 
alf al=Sll ; alf a2=S14 ; alf a4=S44 ; 

'/.Define  elastic  constants  for  constitutive  relations  of  composite-pzt  laminate: 
alf apl=Spll ; alf ap2=Sp!4; alf ap4=Sp44; 
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'/ Define  elastic  constants  for  force-strain  relations  of  composite  laminate  only: 
Bl=l/alf al ; B2=-alf a2/alf al ; B4=alf a4- (alf a2~2/alf al) ; 


'/Define  elastic  constants  for  force-strain  relations  of  composite-pzt  laminate: 
Bpl=l/alf apl;Bp2=-alf ap2/alfapl ;Bp4=alfap4-(alf ap2~2/alf apl) ; 


'/Define  the  elements  of  the  B  matrix: 
bll=4*Bpl*ts*d*T+Bl*t2*h2+Bl*t l*hl ; 
bl2=0 ; 

bl3=-Bl*d*(t2*h2-tl*hl) ; 

b21=0; 

b22=(-Bpl*ts*h2~3*d*T) / (3* (h2-hl) )  +  (1/ (3* (h2-hl) ) ) * (Bpl*ts*hl~3*d*T) - . . . 

(l/12)*Bl*t2*h2~3“(l/12)*Bl*tl*hl~3; 
b23=0 ; 

b31=Bl*d* (t2*h2-tl*hl) ;  b32=0; 

b33= (-4/3) *Bpl*ts*d~3*T-Bl*t2*d~2*h2-Bl*tl*d~2*hl ; 

B=  [bll  bl2  bl3 ;b21  b22  b23;b31  b32  b33] ; 

A=inv(B) ;  all=A(l,l) ;al2=A(l,2) ;al3=A(l,3) ; 
a21=A (2 , 1) ; a22=A (2,2); a23=A (2 , 3) ; 
a31=A (3,1); a32-A (3,2); a33=A (3,3); 

a4=  (1/2)  *d*Bl*t2*h2~2+2*h2*Bpl*ts*d~2*T+d*hl*Bl*t2*h2+4*hl*Bpl*ts*d~2*T+ .  .  . 
(1/2) *d*Bl*tl*hl~2+2*m*d*T*Bl*t2*h2+8*m*Bpl*d~2*T~2*ts ; 

a5=(l/12) *(-Bl*tl*hl~3+6*Bpl*ts*d*T*h2*hl+Bl*t2*h2~3+6*Bpl*ts*d*T*h2~2+ . . . 
8*m*Bpl*d*T~2*ts*h2+16*m*Bpl*d*T's2*ts*hl)  ; 

a6=(-l/2)*d~2*(“4*m*T*Bl*t2*h2“t2*Bpl*h2~2+4*Bl*t2*d*T*h2-4*d*t2*T*Bpl*h2+. . . 
Bl*tl*hl''2-2*Bl*hl*t2*h2)  ; 
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'/.Define  the  elements  of  the  C  matrix: 

cll=B2*h2-B2*hl ;  c21=Bp2*d*T* (h2+hl) ;  c31=B2*d* (h2+hl) ; 

cl2=B2*((-l/2)*Bl*t2*h2~2-2*h2*Bpl*ts*d*T+4*hl*Bpl*ts*d*T+(l/2)*Bl*tl*hl~2+. . . 

hl*Bl*t2*h2)+Bp2* (2*d*T*Bl*t2*h2+4*Bpl*ts*d~2*T~2) ; 
c22=(B2/12)*(Blnl*hl~3+Bl*t2*h2~3)  +  (Bp2/2)*(-d*T*Bl*t2*h2~2-4*Bpl*ts*d~2*,r2*hl-.  .  . 

d*T*hl*Bl*t2*h2-4*Bpl*ts*d~2*T~2*h2) ; 
c32=(-d/6)*(3*B2*Bl*t2*h2~2+12*d*B2*T*Bpl*ts*h2+8*Bpl*ts*d~2*T"2*Bp2+ . . . 
3*B2*Bl*tl*hl~2+6*B2*hl*Bl*t2*h2+24*d*B2*T*hl*Bpl*ts) ; 

cl3=(l/12) *B2* (Bl*t2*h2''3+6*Bpl*ts*h2~2*d*T+6*Bpl*ts*hl*h2*d*T+Bl*tl*hl''3)  ; 
c23=(l/4)*Bpl*ts*Bp2*d~2*T~2*(h2~2+hr2+2*hl*h2)  ; 

c33=(l/12) *d*B2* (-Bl*tl*hl~3+Bl*t2*h2~3+6*Bpl*ts*h2''2*d*T+6*Bpl*ts*hl*h2*d*T) ; 

cl4=(l/6)*d*(-8*Bplns*d~2*’r2*Bp2+3*B2*Bl*tl*hl"2-6*B2*Bl*hl*t2*h2.  .  . 

-12*Bp2*d*T*Bl*t2*h2-12*d*B2*t2*T*h2*Bl+3*B2*t2*Bpl*h2~2+12*d*B2*t2*T*Bpl*h2) ; 
c24= (1/12)  *d*  (B2*Bl*tl*hl~3+6*Bp2*d*T*h2*hl*Bl*t2+6*Bp2*d*T*Bl*t2*h2~2-B2*t2*Bpl*h2~3)  ; 
c34=(l/2)  *d~2*B2*  (-4*d*t2*T*h2*Bl+t2*Bpl*h2~2+4*t2*d*T*Bpl*h2-Bl*tl*hl''2+2*Bl*hl*t2*h2) ; 

C=  [ell  cl2  cl3  cl4; c21  c22  c23  c24;c31  c32  c33  c34] ; 

'/.Define  the  elements  of  the  D  vector: 

dl=(4*Bp4*d*T*tl*t2+B4*ts*tl*h2+B4*ts*t2*hl) / (ts*t2*tl) ; 

d2=(-l/2)  *  (B4*ts*tl*Bl*t2*h2~2+16*Bpl*d'‘2*T''2*Bp4*t2*ts*tl+8*B4*t2*ts~2*hl*Bpl*d*T+ .  .  . 
2*B4*t2~2*ts*Bl*hl*h2+B4*t2*ts*Bl*tl*hl~2+4*Bp4*t2~2*tl*d*T*h2*Bl+ . . . 
4*B4*ts~2*tl*h2*Bpl*d*T) / (t2*ts*tl) ; 

d3= (1/12) * (8*Bpl*d~2*T~2*Bp4n2*h2+16*Bpl*d~2*T~ 2*Bp4*t2*hl+6*B4*d*T*Bpl*ts*h2~2+ . . . 

6*B4*d*T*Bpl*ts*hl*h2+B4*Bl*t2*h2~3-Bl*B4*t2*hl~3) /t2 ; 
d4=(l/2)*d*(4*B4*d*T*ts*tl*Bpl*h2+2*Bp4*Bl*ts*hl*h2*t2-Bp4*Bl*ts*tl*hl~2+. . . 
4*Bp4*tl*d*T*Bl*t2*h2+B4*ts*tl*Bpl*h2~2-4*B4*d*T*ts*tl*h2*Bl)  /  (ts*tl)  ; 


7, Define  the  shear  flow  at  s=0: 

q0= (1/ (2*Area) ) * (M-a22*a5*Vy+ (a!3*a4-a33*a6) *Vx) ; 
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'/.Define  the  derivatives  of  the  strain,  and  curvatures: 
e0p=al3*Vx;  Kxp=a22*Vy;  Kyp=a33*Vx; 

'/.Define  the  values  for  the  Q’s: 

Ql=cll*q0+cl2*e0p+cl3*Kxp+cl4*Kyp; 

Q2=c21*q0+c22*e0p+c23*Kxp+c24*Kyp; 

Q3=c31*q0+c32*e0p+c33*Kxp+c34*Kyp; 

Z1=0 ;  Z2=0 ; Z3=0 ;  eO=al 1* (P-Q1+Z1) -al2* (Mx+Q2+Z2) -al3* (My+Q3+Z3) ; 
Kx=a21* (P-Q1+Z1) -a22* (Mx+Q2+Z2) -a23* (My+Q3+Z3) ; 

Ky=a31* (P-Q1+Z1) -a32* (Mx+Q2+Z2) -a33* (My+Q3+Z3) ; 

'/.The  rotation  in  radians  per  unit  length: 

dFdz=(l/ (2*Area) ) * (-e0*cll+Kx*c21+Ky*c31+q0*dl+e0p*d2+Kxp*d3+Kyp*d4) ; 

'/.The  rotation  in  degrees  per  unit  length: 
dfdz_deg=dFdz* 180/pi ; 

'/.The  rotation  in  degrees  for  entire  length  of  beam: 

Fideg=df dz_deg*L ; 

return 
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Appendix  E.  Anisotropic 


Composite-PZT  Strain  Actuation  Codes 

'/.Program  Compl_aniso_pzt_e  .m 

’/.Written  by  Capt  Peter  Cseke,  Jr.  -  Fall  1999  -  AFIT/ENY 

'/.Driver  program  to  Mat  lab  functions  Layup,  m 

'/,  Anisotorsion_pzt  .m 

’/.Defines  trapezoid  single-cell  box  beam  geometric  and  material  properties, 

’/.applied  cross  section  loads  and  moments. 

’/.Calculates  tip  twist  angle  of  box  beam. 

'/.Uses  :  Layup. m  to  calculate  engineering  properties  of  lamina  and 
*/,  composite-pzt  lamina. 

'/,  Anisotorsion_pzt_e .m  to  calculate  twist  angle  of  box  beam. 

’/.Outputs:  Tip  twist  angle  of  box  beam. 

’/. 

'/.User  enters : 

'/.orientation:  [0  45  -45  90  -90]  etc,  as  a  row  vector; 

'/.times  :  The  number  of  times  ply  is  repeated; 

'/.symmetry  :  0  for  no; 

’/,  1  for  yes 

'/.Example  :  [0  45  -45  90]  4s 
'/,  orient  =  [0  45  -45  90] 

*/,  times  =  4 

*/.  sym  =  1 

'/.Note  :  All  properties  are  vectors  with  properties  per  ply. 

'/.Limits  :  Assumes  weighted  center  of  laminate  is  at  h/2. 

************************************************************  ******  ******** 

close  all;  clear  all; 

’/.Define  beam  geometric  parameters 

E-l 


cr=60.96; 

cl=0.25;c2=0,70; 

hl=8.4334; 

h2=5 . 3035 ; 

121=27.7732; 

122=29.2992; 

L=255 . 96; 


'/.Root  chord  length  [in] 

'/.Main  spar  and  Rear  spar  location  in  chord  percent 
'/.Height  of  Main  Spar  [in] 

'/.Height  of  Rear  Spar  [in] 

'/.Length  of  top  skin  between  Main  and  Rear  Spars  [in] 

7, Length  of  bottom  skin  between  Main  and  Rear  Spars  [in] 
'/.Spanwise  length  of  torque  box  [in] 


’/.Define  the  Loads  and  Moments  acting: 


P=0; 

Vx=0; 

Vy=0; 

Mx=0; 

My=0; 
M_ftlb=0 ; 
M=M_ftlb*12; 


'/.Spanwise  extensional  load  [lb] 
'/.Vertical  shear  load  (lift)  [lb] 
'/.Horizontal  shear  load  (drag)  [lb] 
'/.Applied  moment  about  x  axis  (bending) 
7, Applied  moment  about  y  axis  (bending) 
'/.Applied  moment  about  z  axis  (torque) 

7. 


[in-lb] 

[in-lb] 

[ft-lb] 

[in-lb] 


'/.Voltage  applied  to  PZT  lamina  in  poling  direction 
7.V=[100  250  500  750  1000];  '/.AFC  Lamina 

V=  [50  75  100  150  200];  '/.G-1195  Lamina 


'/.For  the  sake  of  the  legend  for  figure  2,  rewrite  in  string  format 
'/.Uncomment  depending  which  PZT  lamina  is  used 

7,Volts=[’100  V  ’;’250  V  ’  ;’500  V  ’  ;’750  V  ’ ; ’  1000  V1]  ;  '/.AFC  Lamina 
Volts=  [’ 50  V  ’;’75  V  * ;  *  100  V’;’150  V’  1*200  V’]  ;  '/.G-1195  lamina 


</,*****  **  **********************  ****  *******  ************* 
'/.Carbon/Epoxy  (AS4/3501-6) 


El=20.6*10~6; 
E2=l . 50*10~6 ; 
G12=l . 04*10~6; 
vl2=0 . 27 ; 
v21=0 . 02 ; 
thick=0.005; 


'/.Young’s  Modulus  (fiber  direction)  [psi] 
'/.Young’s  Modulus  (transverse  direction)  [psi] 
'/.Shear  Modulus  (in-plane)  [psi] 

'/.Poisson’s  Ratio 
'/.Poisson’s  Ratio 

'/.Thickness  per  layer  of  angle  lamina 
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*/,  (given  by  material  props.) 


y,*  ***********************************************  ****** 
'/.AFC  Piezoceramic  Lamina  (Aaron  Bent) 


'/,Elpzt=4 . 6786*10*6 ; 

'/.E2pzt=2 .4173*10*6; 

’/,G12pzt=5 . 8015*10*5; 
*/.thickpzt_mm=0 .  165 ; 
7,thickpzt_m=thickpzt_mm/ 1000 ; 
*/,thickpzt=thickpzt_m/0 . 0254 ; 
*/,ef  insp__m=0 . 001125 ; 


*/, Young’s  Modulus  (poling,  or  fiber  direction  psi) 
'/.Young’s  Modulus  (transverse  direction)  [psi] 
'/.Shear  Modulus  (in-plane)  [psi] 

'/.Thickness  of  one  layer  of  PZT  lamina  [mm] 

'/.Thickness  of  one  layer  of  PZT  lamina  [m] 

*/. Thickness  of  one  layer  of  PZT  lamina  [in] 

'/.Electrode  finger  spacing  [m] 


'/.thick_select=ef  insp_m; 

'/.d33=180*10*-12; 

'/.d31=-50*10*-12; 


'/.The  AFC  lamina  uses  electrode  f ingerspacing 
*/,as  characteristic  thickness  to  divide  V 
*/, Piezoelectric  constant  in  poling  direction  [m/V] 
'/.Piezoelectric  constant  in  transverse  direction  [m/V] 


*/,PZT-Composite  Lamina  G-1195 
Elpzt=5. 4389*1CT6; 

E2pzt=2 . 0305*10*6; 

G12pzt=5 . 5114*10*5 ; 
thickpzt_mm=0 . 2032 ; 
thickpzt_m=thickpzt_mm/1000; 
thickpzt=thickpzt_m/0 . 0254; 


(Ron  Barrett) 

'/.Young’s  Modulus  (poling  or  fiber  direction  psi) 
'/.Young’s  Modulus  (transverse  direction)  [psi] 
'/.Shear  Modulus  (in-plane,  psi) 

'/.Thickness  of  one  layer  of  PZT  lamina  [mm] 

'/.Thickness  of  one  layer  of  PZT  lamina  [m] 

'/.Thickness  of  one  layer  of  PZT  lamina  [in] 


Lambda=207*10~-6 ;  '/.Actuation  Strain 

thick_select=thickpzt_m;  '/.The  G-1195  lamina  uses  lamina  thickness  to  divide  V'/, 


d33=250*10~-12 ; 
d31=25*10~-12; 


vl2pzt=0 . 30 ; 


'/.Piezoelectric  constant  in  poling  direction  [m/V] 
'/.Piezoelectric  constant  in  transverse  direction  (m/V) 
*/, using  the  small  signal  linear  model 
'/.(Crawley- Anders  on  1990) 

'/.Poisson’s  Ratio 
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anglepzt=45; 


'/.Orientation  of  PZT  patch 


***********♦***************************************************+*********** 
theta=45;  '/.Center  composite  fiber  rotation  angle  (deg) 

times=l;  '/.Number  of  times  lamina  is  wound 

sym=0;  '/.Symmetry  of  lamina  0=no,l=yes 

'/.Define  layer  of  angled  lamina  (only  one  layer  of  variable)  orientation  angles 
anglemin=0 

angleincr=5  '/.Keep  increment  at  5,  so  that  15,  30,  45,  etc.  can  be  found 

anglemax=90; 

'/.The  same  defined  for  easier  plotting 

angle_select=[15  30  45  60  75  90];  anglerange=  [’  15  deg'j'SO 
deg’jMS  deg9;  9  60  deg’;’75  deg’j^O  deg’]  ; 

count=0;  '/.Reset  the  loop-count  variable  (#  of  angles) 

select=l;  '/.Set  counter  for  counting  elements  of  angle-select 

for  j=anglemin:  angle incr :  anglemax  '/.Loop  through  the  angle  range 

count=count+l ;  '/.Update  the  loop  variable 

angle (count) =j ; 

orientals  [0  j  anglepzt  -j  0]; 

'/.Orientation  angles  for  3  layers  of  composite-PZT  laminate: 

*/,orientall=  [0  anglepzt  anglepzt  anglepzt  0]; 

orientcom=[0  j  theta  -j  0];  '/.Orientation  angles  for  composite  laminate 

locatepzt=3;  '/.The  pzt  layer  location  in  the  laminate 

'/.Needs  to  be  changed  to  reflect  changes  in 
'/."orientall"  and  "orientcomp" 

y***************** ******************************* *************************** 
'/.Redefine  thicknesses  according  to  number  of  plies 

t=  [thick];  '/.Thickness  vector  per  layer  of  lamina 
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tl=times*thick*length(orientcom) * (sym+1)  ;  ‘/.Total  main  spar  lamina  thickness 
t2=tl;  '/.Total  rear  spar  lamina  thickness 

*/, Total  thickness  of  laminate  (composite  and  pzt)  at  skin 
ts=times*thick* (length(orientall) -1) * (sym+1) +thickpzt ; 

y* ************************************************************************* 
‘/.Form  the  material  properties  for  each  angle  ply 
‘/.(number  of  angle  ply  from  orientall) 

Elxp=ones (size (orientall) ) *E1 ; 

E2xp=ones (size (orientall) ) *E2 ; 

G12xp=ones (size (orientall) ) *G12; 
v!2xp=ones (size (orientall) ) *vl2 ; 
txp=ones (size (orientall) ) *t ; 

Elxp  (locatepzt )  =Elpzt ;  */,Elxp  (2)  =Elpzt ;  Elxp  (4)  =Elpzt ; 

E2xp  (locatepzt)  =E2pzt ;  '/,E2xp  (2)  =E2pzt ;  E2xp  (4)  =E2pzt  ; 

G12xp  (locatepzt)  =G12pzt ;  */,G12xp  (2)  =G12pzt ;  G12xp(4)  =G12pzt ; 
vl2xp  (locatepzt)  =vl2pzt  ;*/,vl2xp  (2)  =vl2pzt  ;vl2xp  (4)  =vl2pzt ; 
txp  (locatepzt)  =thickpzt ;  */,txp(2)  =thickpzt ;  txp  (4)  =thickpzt ; 

*/,  ************************************************************************** 
‘/.Call  the  Layup. m  subroutine  to  calculate  composite-pzt  laminate  engineering 
‘/.properties 

[Epx,Epy,Gpxy,ap,vpxy,vpyx,npsx,npxs ,npys,npsy]=Layup(Elxp,E2xp, . . . 

vl2xp,G12xp, txp, orientall .times , sym) ; 

*/. The  Compliance  Matrix  (plane  stress)  for  the  Composite-PZT  Laminate  in  the  x-y 
'/.structural  axes : 

Spll=l/Epx;  Spl2=-vpyx/Epy ;  Spl4=npsx/Gpxy ; 

Sp21=-vpxy/Epx ;  Sp22=l/Epy;  Sp24=npsy/Gpxy ; 

Sp41=npxs/Epx ;  Sp42=npys/Epy ;  Sp44=l/Gpxy; 

Spzt= [Spll  Spl2  Spl4; Sp21  Sp22  Sp24;Sp41  Sp42  Sp44] ; 
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'/t  *************************************************************************** 

*/, Calculate  laminate  engineering  properties  (composite  laminate  only) 

Elx=ones (size (orientcom) ) *E1 ; 

E2x=ones (size (orientcom) ) *E2 ; 

G12x=ones (size (orientcom) ) *G12 ; 
vl2x=ones (size (orientcom) ) *vl2 ; 
tx=ones (size (orientcom) ) *t ; 

'/.Call  the  Layup. m  subroutine  to  get  composite  laminate  engineering  properties 
[Ex , Ey , Gxy , a , vxy , vyx , nsx , nxs , nys , nsy] =Layup (Elx , E2x , . . , 

vl2x,G12x,tx, orient com, times , sym) ; 

‘/.The  Compliance  Matrix  (plane  stress)  for  the  Composite  Laminate  in  the 
*/,x-y  structural  axes : 

Sll=l/Ex;  S12=-vyx/Ey;  S14=nsx/Gxy; 

S21=-vxy/Ex ;  S22=l/Ey ;  S24=nsy/Gxy ; 

S41=nxs/Ex;  S42=nys/Ey;  S44=l/Gxy; 

Scom= [Sll  S12  S14; S21  S22  S24;S41  S42  S44] ; 


yt* *************************************************************************** 

for  i=l : length (V) ; 

V (i)  ; 

‘/.The  applied  electric  field  (thick_select  depends  on  type  of  PZT  lamina  used) 
E33=V(i) /thick.select ;  % [V/m] 

y***********3 ft*********  Use  only  for  G-1195  **************************** 

Ea=V(i) /thickpzt_mm  ‘/.Applied  voltage  per  milimeter  [V/mm] 

‘/.From  Crawley  and  Lazarus  1991  using  Figure  Al. 

Lambda=(  (150*lCT-6)/400)  *Ea  7, Induced  Strain  from  linear  approximation 

‘/.From  Crawley  and  Anderson  1990  using  Figure  7. 


E-6 


d33=(  (50*1(T-12)/  (120*10~-6))*  (Lambda) +(240*10~-12) 


'/,d33=Lambda/E33 ;  '/.The  secant  piezoelectric  coefficient  [m/V] 

d31=(l/10)*d33 

********************************+*************************************** 

'/.The  piezoelectric  strain-shear  vector  in  the  material  principle  directions 
pztstrainl2=  [d31  0  0 ;  0  d31  0 ;  0  0  d33]*E33;  */,  {x  y  z}=[T]{l  2  3>  rotation 
*/.pztstrainl2=[d33  0  0 ;  0  d31  0 ;  0  0  d31]*E33;  '/,  {z  x  y}=  [T]  {3  1  2}  rotation 

'/.Subroutine  to  transform  principle  strains/shears  into  structural  strains/shears 
pztstrainzx=transf orm(pztstrainl2 , anglepzt) ; 

'/.The  piezoelectric  strain  and  shear  in  the  z-x  structural  axes 
ep_z=pztstrainzx(3,3) ; 

*/.ep_z=pztstrainzx(l,l)  ; 

gp_zx=pztstrainzx(3, 1) ; 

'/,gp_zx=pztstrainzx(l ,  2) ; 

'/.Call  the  anisotropic  torsion  solution  subroutine  to  calculate  the  amount 
'/.of  tip  torsion  for  the  single  cell  box  beam 

fideg=Anisotorsion_pzt_e(cr,cl,c2,hl,h2,121,122,tl,t2,ts,L,P,Vy,Vx,Mx, . . . 

My , M , ep_z , gp_zx , Spzt , Scorn) ; 

twist  (count  ,i)=fideg;  '/.The  2-D  twist  angle  variable 

'/.Rows:  different  composite  angle 
'/.Columns:  different  Voltage  applied 

end 

'/.Select  only  the  angles  desired  to  plot  (don't  want  to  plot  all  angles  0-90) 
if  j==angle_select (select) 

twist.select (select , : ) =twist (count , : ) ; 
select=select+l ; 
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end 


end 

'/.Here  are  the  angles  from  15-90  with  15  degree  increments 
twist. select 

'/.Declare  plot  style  lines  and  markers 

linestyle= [*b-  *  ;*g-  *  ;*r-  *  ; *m-  * ; ,k-  * ;*c-  *  ; *y-  *]  ; 

markstyle=[*b+  * ; 'go  *  ;*rh  * ; *md  1 ;  *kp  *  ;’c*  * ;*yv  *]  ; 

figure (1) 

for  k=l : 1 : length(angle. select)  '/.Length  of  lamina  angles 

'/.Plot  the  angles  of  twist  with  markers 
plot (V, twist.select (k, : ) ,markstyle(k, : ) ) ; 
hold  on; 

end  for  k^l : 1 : length (angle.select) 

'/.Plot  the  singles  of  twist  with  lines  for  better  visibility 
plot (V, twist .select (k, : ) ,linestyle(k, : ) *) ; 
hold  on; 

end 

'/.title  ( [’Twist  Angles  Due  To  Strain  Actuation  -  Lamina  Angles*]); 

ylabel ( ’Twist  Angle  (deg)*); 

xlabel ( *  Applied  Electric  Field  (V)*); 

legend ( [anglerange] ,2) ; 

grid  on;  hold  off; 

figure (2) 

for  m=l:  1:  length (V)  '/.Length  of  Voltages 

'/.Plot  the  angles  of  twist  with  markers 
plot (angle, twist ( : ,m) ,markstyle (m, : )) ; 
hold  on; 

end  for  m=l : 1 : length (V) 

'/.Plot  the  angles  of  twist  with  lines  for  better  visibility 
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plot (angle , twist ( : ,m) , linestyle (m, : )) ; 
hold  on; 

end 

'/title ( [' Twist  Angles  Due  To  Strain  Actuation  -  Applied  Voltage']); 

ylabel ( 'Twist  Angle  (deg)'); 

xlabel ('Fiber  Angle  (deg)'); 

legend ( [Volts]  ,0) ; 

grid  on;  hold  off; 
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'/.Function  Anisotorsion_pzt_e  .m 


'/.Written  by  Capt  Peter  Cseke,  Jr.  -  Fall  1999  -  AFIT/ENY 

'/.Calculates  the  tip  twist  angle  of  the  single-cell  trapezoid  box-beam  with  PZT  actuator 
'/.embedded  in  top  and  bottom  skins.  No  PZT  in  main  or  rear  spars. 

'/.Uses:  Box  beam  geometric  and  material  parameters  are  inputed  from  driver  program 

'/.Output:  Tip  twist  angle  in  degrees 

'/.Also  works  for  isotropic  torsion,  because  Beta2=0  (B2  and  Bp2)  for  isotropic  materials 

function[Fideg]=Anisotorsion_pzt_e(cr,cl,c2 ,hl ,h2,121 ,122,tl ,t2,ts ,L,P,Vy, Vx,Mx, . . . 

My,M,ep,gp,Sp,Sc) ; 

'/(++:  +  +  *  +  ^  +  +  +  ^*  +  +  *  +  +  +  +  +  +  ++  +  +  ++  +  ***  +  ****+********=f:***************************t* 

'/.Define  the  half -width,  and  top  surface  length  correction  factor 
d=cr*(c2-cl)/2; 

T=sqrt (1+ ( (h2-hl) / (4*d) ) ~2) ; 

Area=2*d* (hl/2+h2/2) ; 
m=(hl+h2)/4; 

'/.Designate  the  compliance  Matrix  for  the  composite  laminate: 

Sll=Sc(l,l) ; S12=Sc (1 ,2) ; S14=Sc (1 , 3) ; 

S21=Sc(2,l) ;S22=Sc(2,2) ;S24=Sc(2,3) ; 

S41=Sc(3,l) ;S42=Sc(3,2) ;S44=Sc(3,3) ; 

'/.Designate  the  compliance  Matrix  for  the  composite-pzt  laminate: 

Spll=Sp(l ,1) ; Spl2=Sp(l ,2) ;Spl4=Sp(l,3) ; 

Sp21=Sp (2,1) ; Sp22=Sp (2 , 2) ; Sp24=Sp (2 , 3) ; 

Sp41=Sp(3,l) ;Sp42-Sp(3,2) ;Sp44=Sp(3,3) ; 

'/.Define  elastic  constants  for  constitutive  relations  of  composite  laminate  only: 
alf al=Sll ; alf a2=S14 ; alf a4=S44 ; 

*/, Define  elastic  constants  for  constitutive  relations  of  composite-pzt  laminate: 
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alf apl=Spll ; alf ap2=Spl4 ; alf ap4=Sp44; 


'/♦Define  elastic  constants  for  force-strain  relations  of  composite  laminate  only 
Bl=l/alfal;B2=-alfa2/alfal;B4=alfa4-(alfa2~2/alfal) ; 

'/♦Define  elastic  constants  for  force-strain  relations  of  composite-pzt  laminate: 
Bpl=l/ alfapl ;Bp2=-alf ap2/alf apl ;Bp4=alf ap4- (alf ap2~2/ alf apl) ; 

'/♦Define  the  elements  of  the  B  matrix: 
bll=4*Bpl*ts*d*T+Bl*t2*h2+Bl*tl*hl ;  bl2=0; 
bl3=-Bl*d* (t2*h2-tl*hl) ; 


b21=0; 

b22=(-Bpl*ts*h2~3*d*T) / (3* (h2-hl) ) + (1/ (3* (h2-hl) ) ) * (Bpl*ts*hl~3*d*T) - . . . 

(1/12) *Bl*t2*h2~3- (1/12) *Bl*tl*hl~3 ; 
b23=0; 

b31=Bl*d*(t2*h2-tl*hl) ; 
b32=0 ; 

b33=(-4/3) *Bpl*ts*d~3*T-Bl*t2*d~2*h2-Bl*tl*d~2*hl; 

B=  [bll  bl2  bl3;b21  b22  b23;b31  b32  b33] ; 

A=inv(B) ; 

all=A(l , 1) ; al2=A(l ,2) ; al3=A(l ,3) ; 
a21=A(2, 1) ; a22=A(2,2) ;a23=A(2,3) ; 
a31=A(3, 1) ; a32=A(3 , 2) ; a33=A(3,3) ; 

a4=(l/2)  *d*Bl*t2*h2~2+2*h2*Bpl*ts*d'N2*T+d*hl*Bl*t2*h2+4*hl*Bpl*ts*d~2*T+ .  .  . 
(1/2) *d*Bl*tl*hl~2+2*m*d*T*Bl*t2*h2+8*m*Bpl*d~2*T~2*ts ; 

a5=(l/12)  *  (-Bl*tl*hl~3+6*Bpl*ts*d*T*h2*hl+Bl*t2*h2~3+6*Bpl*ts*d*T*h2~2+ .  .  . 
8*m*Bpl*d*T~2*ts*h2+16*m*Bpl*d*T~2*ts*hl) ; 
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a6=(-l/2)*d~2*(-4*m*T*Bl*t2*h2-t2*Bpl*h2~2+4*Bl*t2*d*T*h2-4*d*t2*T*Bpl*h2+. . . 
Bl*tl*hl"2-2*Bl*hl*t2*h2) ; 

*/, Define  the  elements  of  the  C  matrix: 
cll=B2*h2-B2*hl ; 
c21=Bp2*d*T* (h2+hl) ; 
c31=B2*d* (h2+hl) ; 

cl2=B2* ( (-1/2) *Bln2*h2~2-2*h2*Bpl*ts*d*T+4*hl*Bpl*ts*d*T+  (1/2) *Bl*t l*hl~2+ .  .  . 

hl*Bl*t2*h2)+Bp2*  (2*d*T^*t2*h2+4*Bpl*ts*d~2*T''2) ; 
c22=  (B2/12)  *  (Bl*tl*hl~3+Bl*t2*h2~3)  +  (Bp2/2)  *  (-d*T*Bl*t2*h2~2-4*Bpl*ts*d'‘2*T~2*hl- .  .  . 

d*T*hl*Bl*t2*h2-4*Bpl*ts*d~2*T~2*h2) ; 
c32=(-d/6)*(3*B2*Bl*t2*h2~2+12*d*B2*T*Bpl*ts*h2+8*Bpl*ts*d~2*T"2*Bp2+. . . 
3*B2*Bl*tl*hl''2+6*B2*hl*Bl*t2*h2+24*d*B2*T*hl*Bpl*ts)  ; 

cl3=(l/12)*B2*(Bl*t2*h2~3+6*Bpms*h2~2*d*T+6*Bpl*ts*hl*h2*d*T+Bl*tl*hl~3)  ; 
c23=(l/4) *Bpl*ts*Bp2*d~2*T~2*(h2~2+hl~2+2*hl*h2) ; 

c33=(l/12)  *d*B2*  (-Bl*tl*hl~3+Bl*t2*h2~3+6*Bpl*ts*h2~2*d*T+6*Bpl*ts*hl*h2*d*T) ; 

c  14=  (1/6)  *d*  (-8*Bpl*ts*d''2*T~2*Bp2+3*B2*Bl*tl*hl~2-6*B2*Bl*hl*t2*h2 .  .  . 

-12*Bp2*d*T*Bl*t2*h2-12*d*B2*t2*T*h2*Bl+3*B2*t2*Bpl+h2~2+12+d*B2*t2*T*Bpl*h2) ; 
c24=(l/12)*d*(B2*Bl*tl*hl~3+6*Bp2*d*T*h2*hl*Bl*t2+6*Bp2*d*T*Bl*t2*h2~2-B2*t2*Bpl*h2~3) ; 
c34=(l/2)  *d"2*B2*  (-4*d*t2*T*h2*Bl+t2*Bpl*h2'‘2+4*t2*d*T*Bpl*h2-Bl*tl*hl~2+2*Bl*hl*t2*h2)  ; 

C=  [ell  cl2  c!3  cl4 ; c21  c22  c23  c24;c31  c32  c33  c34] ; 


*/, Define  the  elements  of  the  D  vector: 
dl=(4*Bp4*d*T*tl*t2+B4*ts*tl*h2+B4*ts*t2*hl) / (ts*t2*tl) ; 

d2=(-l/2)  *  (B4*ts*tl*Bl*t2*h2~2+16*Bpl*d~2*T~2*Bp4*t2*ts*tl+8*B4*t2*ts',2*hl*Bpl*d*T+ .  .  . 
2*B4*t2~2*ts*Bl*hl*h2+B4*t2*ts*Bl*tl*hl''2+4*Bp4*t2~2*tl*d*T*h2*Bl+ .  .  . 
4*B4*ts~2*tl*h2*Bpl*d*T)/(t2*ts*tl) ; 

d3=(l/12)*(8*Bpl*d~2*T~2*Bp4*t2*h2+16*Bpl*d~2*T~2*Bp4*t2*hl+6*B4*d*T*Bpl*ts*h2~2+. . . 

6*B4*d*T*Bpl*ts*hl*h2+B4*Bl*t2*h2~3-Bl*B4*t2*hl~3) /t2 ; 
d4=(l/2) *d* (4*B4*d*T*ts*tl*Bpl*h2+2*Bp4*Bl*ts*hl*h2*t2-Bp4*Bl*ts*tl*hl~2+ . . . 
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4*Bp4*tl*d*T*Bl*t2*h2+B4ns*tl*Bpl*h2~2-4*B4*d*T*ts*tl*h2*Bl)/(ts*tl)  ; 


'/, Define  the  shear  flow  at  s=0: 

qO=(l/ (2*Area) ) * (M-a22*a5*Vy+(al3*a4-a33*a6) *Vx) ; 

'/, Define  the  derivatives  of  the  strain,  and  curvatures: 
e0p=al3*Vx; 

Kxp=a22*Vy ; 

Kyp=a33*Vx; 

*/, Define  the  values  for  the  QJs: 
Ql=cll*q0+cl2*e0p+cl3*Kxp+cl4*Kyp; 
Q2=c21*q0+c22*e0p+c23*Kxp+c24*Kyp; 
Q3=c31*q0+c32*e0p+c33*Kxp+c34*Kyp; 


Zl=4*ep*d*T*Bpl*ts ; 

Z2=0;  Z3=0; 

eO=all* (P-Q1+Z1) -al2* (Mx+Q2-Z2) -al3* (My+Q3-Z3) ; 

Kx=a21* (P-Q1+Z1) -a22* (Mx+Q2-Z2) -a23* (My+Q3-Z3) ; 

Ky=a31* (P-Q1+Z1) -a32* (Mx+Q2-Z2) -a33* (My+Q3-Z3) ; 

’/.The  PZT  shear  around  the  perimeter  (using  shear  in  the  structural  axes) 
g_total=4*gp*d*T ; 
e_total=4*ep*d*T ; 

’/.The  rotation  in  radians  per  unit  length: 

dFdz=(l/(2*Area) ) * (-e0*cll+Kx*c21+Ky*c31+q0*dl+e0p*d2+Kxp*d3+Kyp*d4+e_total*cll+g_total) ; 

*/,The  rotation  in  degrees  per  unit  length: 
df dz_deg=dFdz*180/pi ; 

'/.The  rotation  in  degrees  for  entire  length  of  beam: 

Fideg=df dz_deg*L ; 
return 
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■/.Function  transform.m 


7, Written  by  Capt  Peter  Cseke,  Jr.  -  Winter  2000  -  AFIT/ENY 
/.Subroutine  to  Compl_aniso_pzt_e  .m 

X 


'/.Calcultes  the  transformed  values  of  the  piezoelectric  strain  vector  in  the 
*/, structural  axes,  given  the  PZT  layup  angle. 


7. 

'/.Caution: 

7. 

7. 

X 

X 

X 

X 

X 

X 

'/.Uses  : 
'/Output  s : 


Fiber  orientation  (in  this  subroutine  only)  agrees  with  that 
accepted  in  the  literature  for  smart  materials,  that  is: 

3  -  fiber  direction 

1  -  transverse  direction 

2  -  out  of  plane  direction 

This  change  in  notation  does  not  affect  the  notation  in  main  code 
as  long  as  the  PZT  strains  are  sent  back  in  the  z-x  structural  axes 
that  correspond  to  the  original  structural  axes  defined  in  the 
problem  statement . 

PZT  strains  (3x3  matrix)  in  principle  axes  1-2 
PZT  strains  (3x3  matrix)  in  structural  axes  z-x. 


function  [pztstrainzx] =transf orm(pztstrain!2 , angledeg) ; 


'/.Calculate  the  rotation  angle  in  radians 
anglerad=angledeg*pi/180 ; 


'/.For  the  sake  of  simplicity 
s=s in (angler ad) ; 
c=cos (anglerad) ; 


'/.The  transformation  matrix 

Trans=  [c  0  s ;  0  1  0 ; -s  0  c]  ;  '/,  {x  y  z}=[T]{l 
*/,Trans=[c  -s  0;s  c  0;0  0  1];  7.  {z  x  y}=[T]{3 


2  3}  rotation 
1  2}  rotation 
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'/.Invert  the  transformation  matrix 
Cinv=inv (Trans) ; 

'/.The  transformation  of  PZT  strains  from  the  3-1  principle  into  the  z-x  structural  axes 
pztstrainzx=Cinv*pztstrainl2*Trans ; 

return 
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